diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index ba02a82..fdcb201 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -7,42 +7,125 @@ namespace polyMPO{ void printVTP_mesh(MPMesh& mpMesh, int printVTPIndex=-1); -void MPMesh::calcBasis() { - assert(p_mesh->getGeomType() == geom_spherical_surf); - auto MPsPosition = p_MPs->getPositions(); - auto mp_basis_field = p_MPs->getData(); // we can implement MPs->getBasisVals() like MPs->getPositions() - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - auto vtxCoords = p_mesh->getMeshField(); - double radius = p_mesh->getSphereRadius(); - - auto calcbasis = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { - if(mask) { //if material point is 'active'/'enabled' - Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2)); - // formating - Vec3d v3d[maxVtxsPerElm+1]; - int numVtx = elm2VtxConn(elm,0); - for(int i = 1; i<=numVtx; i++){ - v3d[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); - v3d[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); - v3d[i-1][2] = vtxCoords(elm2VtxConn(elm,i)-1,2); - } - v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); - v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); - v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2); - - double basisByArea3d[maxVtxsPerElm] = {0.0}; - initArray(basisByArea3d,maxVtxsPerElm,0.0); +void MPMesh::calculateStrain(){ + auto MPsPosition = p_MPs->getPositions(); + auto MPsBasis = p_MPs->getData(); + auto MPsAppID = p_MPs->getData(); + + //Mesh Fields + auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField(); + auto gnomProjVtx = p_mesh->getMeshField(); + auto gnomProjElmCenter = p_mesh->getMeshField(); + auto elm2VtxConn = p_mesh->getElm2VtxConn(); + auto velField = p_mesh->getMeshField(); + double radius = 1.0; + if(p_mesh->getGeomType() == geom_spherical_surf) + radius=p_mesh->getSphereRadius(); + + bool isRotated = p_mesh->getRotatedFlag(); + + auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ + if(mask){ + int numVtx = elm2VtxConn(elm,0); + + Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2)); + if(isRotated){ + position3d[0] = -MPsPosition(mp, 2); + position3d[2] = MPsPosition(mp, 0); + } + auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL); + double mpProjX, mpProjY; + computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY); + auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL); + + double basisByArea[maxVtxsPerElm] = {0.0}; + initArray(basisByArea,maxVtxsPerElm,0.0); + double gradBasisByArea[2*maxVtxsPerElm] = {0.0}; + initArray(gradBasisByArea,maxVtxsPerElm,0.0); + + wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, radius, basisByArea, gradBasisByArea); + + double v11 = 0.0; + double v12 = 0.0; + double v21 = 0.0; + double v22 = 0.0; + double uTanOverR = 0.0; + double vTanOverR = 0.0; + for (int i = 0; i < numVtx; i++){ + int iVertex = elm2VtxConn(elm, i+1)-1; + v11 = v11 + gradBasisByArea[i*2 + 0] * velField(iVertex, 0); + v12 = v12 + gradBasisByArea[i*2 + 1] * velField(iVertex, 0); + v21 = v21 + gradBasisByArea[i*2 + 0] * velField(iVertex, 1); + v22 = v22 + gradBasisByArea[i*2 + 1] * velField(iVertex, 1); + uTanOverR = uTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 0); + vTanOverR = vTanOverR + basisByArea[i] * tanLatVertexRotatedOverRadius(iVertex, 0) * velField(iVertex, 1); + } + } + }; + p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate"); +} - // calc basis - getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d); - - // fill step - for(int i=0; i<= numVtx; i++){ - mp_basis_field(mp,i) = basisByArea3d[i]; - } - } - }; - p_MPs->parallel_for(calcbasis, "calcbasis"); +void MPMesh::calcBasis() { + assert(p_mesh->getGeomType() == geom_spherical_surf); + + auto MPsPosition = p_MPs->getPositions(); + auto MPsBasis = p_MPs->getData(); + auto MPsAppID = p_MPs->getData(); + + auto elm2VtxConn = p_mesh->getElm2VtxConn(); + auto vtxCoords = p_mesh->getMeshField(); + double radius = 1.0; + if(p_mesh->getGeomType() == geom_spherical_surf) + radius=p_mesh->getSphereRadius(); + + //For Gnomonic Projection + auto gnomProjVtx = p_mesh->getMeshField(); + auto gnomProjElmCenter = p_mesh->getMeshField(); + + bool isRotated = p_mesh->getRotatedFlag(); + + auto calcbasis = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { + if(mask) { //if material point is 'active'/'enabled' + int numVtx = elm2VtxConn(elm,0); + Vec3d position3d(MPsPosition(mp, 0),MPsPosition(mp, 1),MPsPosition(mp, 2)); + if(isRotated){ + position3d[0] = -MPsPosition(mp, 2); + position3d[2] = MPsPosition(mp, 0); + } + + double mpProjX, mpProjY; + auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL); + computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY); + auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL); + + double basisByArea[maxVtxsPerElm] = {0.0}; + initArray(basisByArea,maxVtxsPerElm, 0.0); + double gradBasisByArea[2*maxVtxsPerElm] = {0.0}; + initArray(gradBasisByArea,maxVtxsPerElm, 0.0); + + wachpress_weights_grads_2D(numVtx, gnom_vtx_subview, mpProjX, mpProjY, radius, basisByArea, gradBasisByArea); + + for(int i=0; i<= numVtx; i++){ + MPsBasis(mp,i) = basisByArea[i]; + } + + //Old method where basis functions calculated using 3D Area + /* + Vec3d v3d[maxVtxsPerElm+1]; + int numVtx = elm2VtxConn(elm,0); + for(int i = 1; i<=numVtx; i++){ + v3d[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); + v3d[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); + v3d[i-1][2] = vtxCoords(elm2VtxConn(elm,i)-1,2); + } + v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); + v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); + v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2); + getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea); + */ + } + }; + p_MPs->parallel_for(calcbasis, "calcbasis"); } void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){ @@ -51,10 +134,10 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){ auto elm2VtxConn = p_mesh->getElm2VtxConn(); auto elm2ElmConn = p_mesh->getElm2ElmConn(); auto MPs2Elm = p_MPs->getData(); - const auto vtxCoords = p_mesh->getMeshField(); + const auto vtxCoords = p_mesh->getMeshField(); auto mpPositions = p_MPs->getData(); Kokkos::View edgeCenters("EdgeCenters",numElms); - + Kokkos::parallel_for("calcEdgeCenter", numElms, KOKKOS_LAMBDA(const int elm){ int numVtx = elm2VtxConn(elm,0); int v[maxVtxsPerElm]; @@ -67,7 +150,7 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){ edgeCenters(elm,i) = (v_ip1 + v_i)*0.5; } }); - + auto CVTEdgeTracking = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ Vec2d MP(mpPositions(mp,0),mpPositions(mp,1));//XXX:the input is XYZ, but we only support 2d vector if(mask){ @@ -90,7 +173,7 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){ } if(currentDistSq < minDistSq){ edgeIndex = i+1; - minDistSq = currentDistSq; + minDistSq = currentDistSq; } } if(edgeIndex <0){ @@ -149,13 +232,14 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ Vec3d center(elmCenter(iElm, 0), elmCenter(iElm, 1), elmCenter(iElm, 2)); Vec3d delta = MPnew - center; - + double minDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2]; int closestElm = -1; //go through all the connected elm, calc distance for(int i=1; i<=numConnElms; i++){ int elmID = elm2ElmConn(iElm,i)-1; - + if (elmID >= numElms) + continue; //New delta Vec3d center(elmCenter(elmID, 0), elmCenter(elmID, 1), elmCenter(elmID, 2)); delta = MPnew - center; @@ -243,7 +327,7 @@ void MPMesh::T2LTracking(Vec2dView dx){ auto mpPositions = p_MPs->getData(); auto MPs2Elm = p_MPs->getData(); auto mpStatus = p_MPs->getData(); - + auto T2LCalc = PS_LAMBDA(const int& elm, const int& mp, const int&mask){ Vec2d MP(mpPositions(mp,0),mpPositions(mp,1));//XXX:the input is XYZ, but we only support 2d vector if(mask){ @@ -299,15 +383,13 @@ void MPMesh::T2LTracking(Vec2dView dx){ } void MPMesh::reconstructSlices() { - if (reconstructSlice.size() == 0) return; - Kokkos::Timer timer; - calcBasis(); - resetPreComputeFlag(); - for (auto const& [index, reconstruct] : reconstructSlice) { - if (reconstruct) reconstruct(); - } - reconstructSlice.clear(); - pumipic::RecordTime("PolyMPO_Reconstruct", timer.seconds()); + if (reconstructSlice.size() == 0) return; + Kokkos::Timer timer; + for (auto const& [index, reconstruct] : reconstructSlice) { + if (reconstruct) reconstruct(); + } + reconstructSlice.clear(); + pumipic::RecordTime("PolyMPO_Reconstruct", timer.seconds()); } bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) { @@ -329,12 +411,13 @@ void MPMesh::push_ahead(){ //Latitude Longitude increment at mesh vertices and interpolate to particle position p_mesh->computeRotLatLonIncr(); - //Interpolates latitude longitude increments and mesh velocity increments to - //MP positions - sphericalInterpolationDispVelIncr(*this); + //Interpolates latitude longitude, mesh velocity increments to MPs + calcBasis(); + sphericalInterpolation(*this); + sphericalInterpolation(*this); //Push the MPs - p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); + p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius(), p_mesh->getRotatedFlag()); pumipic::RecordTime("PolyMPO_interpolateAndPush", timer.seconds()); } @@ -361,16 +444,14 @@ void MPMesh::push_swap_pos(){ p_MPs->updateMPSlice(); } - -void MPMesh::push(){ - +void MPMesh::push(){ Kokkos::Timer timer; p_mesh->computeRotLatLonIncr(); sphericalInterpolation(*this); - p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ + p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius(), p_mesh->getRotatedFlag()); auto elm2Process = p_mesh->getElm2Process(); @@ -396,43 +477,43 @@ void MPMesh::push(){ } void MPMesh::printVTP_mesh(int printVTPIndex){ - auto vtxCoords = p_mesh->getMeshField(); - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - - auto MPsPosition = p_MPs->getPositions(); - - char* fileOutput = (char *)malloc(sizeof(char) * 256); - sprintf(fileOutput,"polyMPO_MPMesh_mesh_%d.vtp", printVTPIndex); - FILE * pFile = fopen(fileOutput,"w"); - free(fileOutput); - - auto h_vtxCoords = Kokkos::create_mirror_view(vtxCoords); - IntVtx2ElmView::HostMirror h_elm2VtxConn = Kokkos::create_mirror_view(elm2VtxConn); - const int nCells = p_mesh->getNumElements(); - const int nVertices = p_mesh->getNumVertices(); - Kokkos::deep_copy(h_vtxCoords,vtxCoords); - Kokkos::deep_copy(h_elm2VtxConn,elm2VtxConn); - fprintf(pFile, "\n \n \n \n \n",nVertices,nCells); - for(int i=0; i\n \n \n \n"); - for(int i=0; i\n \n"); + auto vtxCoords = p_mesh->getMeshField(); + auto elm2VtxConn = p_mesh->getElm2VtxConn(); + + auto MPsPosition = p_MPs->getPositions(); + + char* fileOutput = (char *)malloc(sizeof(char) * 256); + sprintf(fileOutput,"polyMPO_MPMesh_mesh_%d.vtp", printVTPIndex); + FILE * pFile = fopen(fileOutput,"w"); + free(fileOutput); + + auto h_vtxCoords = Kokkos::create_mirror_view(vtxCoords); + IntVtx2ElmView::HostMirror h_elm2VtxConn = Kokkos::create_mirror_view(elm2VtxConn); + const int nCells = p_mesh->getNumElements(); + const int nVertices = p_mesh->getNumVertices(); + Kokkos::deep_copy(h_vtxCoords,vtxCoords); + Kokkos::deep_copy(h_elm2VtxConn,elm2VtxConn); + fprintf(pFile, "\n \n \n \n \n",nVertices,nCells); + for(int i=0; i\n \n \n \n"); + for(int i=0; i\n \n"); - int count = 0; - for(int i=0;i\n \n \n \n\n"); - fclose(pFile); + int count = 0; + for(int i=0;i\n \n \n \n\n"); + fclose(pFile); } } diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp index cb5c893..798c9e2 100644 --- a/src/pmpo_MPMesh.hpp +++ b/src/pmpo_MPMesh.hpp @@ -18,30 +18,36 @@ template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_OnSurfVeloI #define maxMPsPerElm 8 class MPMesh{ - private: - - bool isPreComputed; - + public: - - MPMesh() : isPreComputed(false){}; - void computeMatricesAndSolve(); - void resetPreComputeFlag(); - Kokkos::View precomputedVtxCoeffs; Mesh* p_mesh; MaterialPoints* p_MPs; - std::map> reconstructSlice = std::map>(); - + //For MPI Communication + int numOwnersTot, numHalosTot; + std::vector numOwnersOnOtherProcs; + std::vector numHalosOnOtherProcs; + std::vectorhaloOwnerProcs; + std::vector> haloOwnerLocalIDs; + std::vector> ownerOwnerLocalIDs; + std::vector> ownerHaloLocalIDs; + + void startCommunication(); + void communicateFields(const std::vector>& fieldData, const int numEntities, const int numEntries, int mode, + std::vector>& recvIDVec, std::vector>& recvDataVec); + void communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, int numEntries, int mode, int op); + MPMesh(Mesh* inMesh, MaterialPoints* inMPs): - p_mesh(inMesh), p_MPs(inMPs) { + p_mesh(inMesh), p_MPs(inMPs) { }; + ~MPMesh() { delete p_mesh; delete p_MPs; } + //MP advection and tracking void CVTTrackingEdgeCenterBased(Vec2dView dx); void CVTTrackingElmCenterBased(const int printVTPIndex = -1); void T2LTracking(Vec2dView dx); @@ -50,38 +56,37 @@ class MPMesh{ void push_swap(); void push_swap_pos(); void push(); + + //Used before advection to interpolate fields from mesh to MPs + //And also before reconstruction void calcBasis(); + //Reconstruction DoubleView assemblyV0(); - template - DoubleView wtScaAssembly(); - template - Vec2dView wtVec2Assembly(); - template - void assembly(int order, MeshFieldType type, bool basisWeightFlag, bool massWeightFlag); template void assemblyVtx0(); template void assemblyElm0(); template void assemblyVtx1(); - template - void subAssemblyVtx1(int vtxPerElm, int nCells, int comp, double* array); - - void subAssemblyCoeffs(int vtxPerElm, int nCells, double* m11, double* m12, double* m13, double* m14, - double* m22, double* m23, double* m24, - double* m33, double* m34, - double* m44); - void solveMatrixAndRegularize(int nVertices, double* m11, double* m12, double* m13, double* m14, - double* m22, double* m23, double* m24, - double* m33, double* m34, - double* m44); + void reconstruct_coeff_full(); + void invertMatrix(const Kokkos::View& vtxMatrices, const double& radius); + Kokkos::View precomputedVtxCoeffs_new; + //Not used currently + std::map> reconstructSlice = std::map>(); + template + DoubleView wtScaAssembly(); + template + Vec2dView wtVec2Assembly(); + template + void assembly(int order, MeshFieldType type, bool basisWeightFlag, bool massWeightFlag); template void setReconstructSlice(int order, MeshFieldType type); void reconstructSlices(); void printVTP_mesh(int printVTPIndex); + void calculateStrain(); }; }//namespace polyMPO end diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp index 033386e..2a9af6e 100644 --- a/src/pmpo_MPMesh_assembly.hpp +++ b/src/pmpo_MPMesh_assembly.hpp @@ -6,34 +6,31 @@ namespace polyMPO{ DoubleView MPMesh::assemblyV0(){ - int numVtxs = p_mesh->getNumVertices(); - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - - DoubleView vField("vField2",numVtxs); - auto mpPositions = p_MPs->getData(); //get the array of MP coordinates/positions - auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { - //for elm in elementsInMesh { //pseudo code - the 'parallel_for' handles this - // for mp in materialPointsInElm { //pseudo code (cont.) - if(mask) { //if material point is 'active'/'enabled' - int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element - for(int i=0; iparallel_for(assemble, "assembly"); - return vField; + int numVtxs = p_mesh->getNumVertices(); + auto elm2VtxConn = p_mesh->getElm2VtxConn(); + + DoubleView vField("vField2",numVtxs); + auto mpPositions = p_MPs->getData(); //get the array of MP coordinates/positions + auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { + if(mask) { //if material point is 'active'/'enabled' + int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element + for(int i=0; iparallel_for(assemble, "assembly"); + return vField; } template void MPMesh::assemblyVtx0(){ Kokkos::Timer timer; + constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice; - auto elm2VtxConn = p_mesh->getElm2VtxConn(); + auto elm2VtxConn = p_mesh->getElm2VtxConn(); auto mpData = p_MPs->getData(); const int numEntries = mpSliceToNumEntries(); @@ -91,44 +88,45 @@ void MPMesh::assemblyElm0() { Kokkos::MDRangePolicy> policy({0,0},{numElms, numEntries}); Kokkos::parallel_for("assembly average", policy, KOKKOS_LAMBDA(const int elm, const int entry){ - if (mpsPerElm(elm) > 0){ + if (mpsPerElm(elm) > 0){ meshField(elm, entry) /= mpsPerElm(elm); - } - }); + } + }); pumipic::RecordTime("PolyMPO_Reconstruct_Elm0", timer.seconds()); } - -void MPMesh::resetPreComputeFlag(){ - isPreComputed = false; -} - -//Method 1 for coefficients -void MPMesh::computeMatricesAndSolve(){ +void MPMesh::reconstruct_coeff_full(){ + std::cout<<__FUNCTION__<getMPIComm(); + MPI_Comm_rank(comm, &self); + MPI_Comm_size(comm, &numProcsTot); + //Mesh Information - auto elm2VtxConn = p_mesh->getElm2VtxConn(); + auto elm2VtxConn = p_mesh->getElm2VtxConn(); int numVtx = p_mesh->getNumVertices(); auto vtxCoords = p_mesh->getMeshField(); - + int numVertices = p_mesh->getNumVertices(); //Dual Element Area for Regularization auto dual_triangle_area=p_mesh->getMeshField(); //Material Points + calcBasis(); + auto weight = p_MPs->getData(); - auto mpPositions = p_MPs->getData(); + auto mpPos = p_MPs->getData(); //Matrix for each vertex - Kokkos::View VtxMatrices("VtxMatrices", p_mesh->getNumVertices()); + constexpr int numEntriesMatrix=10; + Kokkos::View vtxMatrices("VtxMatrices", p_mesh->getNumVertices()); + Kokkos::deep_copy(vtxMatrices, 0); //Earth Radius double radius = 1.0; if(p_mesh->getGeomType() == geom_spherical_surf) radius=p_mesh->getSphereRadius(); - bool scaling=true; - int reg_method = 2; - //Assemble matrix for each vertex auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { if(mask) { //if material point is 'active'/'enabled' @@ -136,316 +134,170 @@ void MPMesh::computeMatricesAndSolve(){ for(int i=0; iparallel_for(assemble, "assembly"); - - //Assemble matrix for each vertex and apply regularization - Kokkos::View VtxCoeffs("VtxCoeffs", p_mesh->getNumVertices()); - - Kokkos::parallel_for("solving Ax=b", numVtx, KOKKOS_LAMBDA(const int vtx){ - Vec4d v0 = {VtxMatrices(vtx,0,0), VtxMatrices(vtx,0,1), VtxMatrices(vtx,0,2), VtxMatrices(vtx,0,3)}; - Vec4d v1 = {VtxMatrices(vtx,1,0), VtxMatrices(vtx,1,1), VtxMatrices(vtx,1,2), VtxMatrices(vtx,1,3)}; - Vec4d v2 = {VtxMatrices(vtx,2,0), VtxMatrices(vtx,2,1), VtxMatrices(vtx,2,2), VtxMatrices(vtx,2,3)}; - Vec4d v3 = {VtxMatrices(vtx,3,0), VtxMatrices(vtx,3,1), VtxMatrices(vtx,3,2), VtxMatrices(vtx,3,3)}; - //Define the matrices - Matrix4d A = {v0,v1,v2,v3}; - Matrix4d A_regularized = {v0, v1, v2, v3}; - //Regularization - switch(reg_method){ - case 0:{ - break; - } - case 1:{ - double A_trace = A.trace(); - A_regularized.addToDiag(A_trace*1e-8); - break; - } - case 2:{ - double regParam=sqrt(EPSILON)*(VtxMatrices(vtx,0,0)+VtxMatrices(vtx,1,1)+ - VtxMatrices(vtx,2,2)+VtxMatrices(vtx,3,3)); - A_regularized.addToDiag(regParam); - break; - } - default:{ - printf("Invalid regularization method \n"); - break; - } - } - //Solve Ax=b - double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0}; - CholeskySolve4d_UnitRHS(A_regularized, coeff); - // Undo scaling - double mScale=1; - if(scaling) - mScale=sqrt(dual_triangle_area(vtx,0))/radius; - - coeff[0]=coeff[0]*mScale*mScale; - coeff[1]=coeff[1]*mScale; - coeff[2]=coeff[2]*mScale; - coeff[3]=coeff[3]*mScale; - for (int i=0; iprecomputedVtxCoeffs = VtxCoeffs; - pumipic::RecordTime("PolyMPO_Calculate_MLS_Coeff", timer.seconds()); -} + Kokkos::fence(); + pumipic::RecordTime("Assemble Matrix Per Process" + std::to_string(self), timer.seconds()); + //Mode 0 is Gather: Halos Send to Owners + //Mode 1 is Scatter: Owners Send to Halos + //Op 0 is addition + //Op 1 is replacement + timer.reset(); + int mode = 0; + int op = 0; + if (numProcsTot >1){ + communicate_and_take_halo_contributions(vtxMatrices, numVertices, numEntriesMatrix, mode, op); + mode=1; + op=1; + communicate_and_take_halo_contributions(vtxMatrices, numVertices, numEntriesMatrix, mode, op); + } + pumipic::RecordTime("Communicate Matrix Values" + std::to_string(self), timer.seconds()); -//Method 2 for coefficients -void MPMesh::subAssemblyCoeffs(int vtxPerElm, int nCells, double* m11, double* m12, double* m13, double* m14, - double* m22, double* m23, double* m24, - double* m33, double* m34, - double* m44){ - - Kokkos::Timer timer; - //Material Points Information - MPI_Comm comm = p_MPs->getMPIComm(); - int comm_rank; - MPI_Comm_rank(comm, &comm_rank); + invertMatrix(vtxMatrices, radius); +} - //Mesh Information - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - int numVtx = p_mesh->getNumVertices(); +void MPMesh::invertMatrix(const Kokkos::View& vtxMatrices, const double& radius){ + std::cout<<__FUNCTION__<getNumVertices(); auto vtxCoords = p_mesh->getMeshField(); - auto elm2Process = p_mesh->getElm2Process(); - auto elm2global = p_mesh->getElmGlobal(); + auto dual_triangle_area = p_mesh->getMeshField(); + bool isRotated = p_mesh->getRotatedFlag(); - //Dual Element Area for Regularization - auto dual_triangle_area=p_mesh->getMeshField(); - - //Material Points - calcBasis(); - auto weight = p_MPs->getData(); - auto mpPos = p_MPs->getData(); - auto mpAppID = p_MPs->getData(); + double eps = 1e-7; + double truncateFactor = 0.05; - //Radius - double radius = 1.0; - if(p_mesh->getGeomType() == geom_spherical_surf) - radius=p_mesh->getSphereRadius(); + Kokkos::View VtxCoeffs("VtxCoeffs", nVertices); + Kokkos::deep_copy(VtxCoeffs, 0.0); + + Kokkos::parallel_for("invertMatrix", nVertices, KOKKOS_LAMBDA(const int vtx){ + if(vtxMatrices(vtx, 0) < eps) + return; + + auto small = eps * vtxMatrices(vtx, 0) * dual_triangle_area(vtx, 0)/(radius*radius); + auto truncate = truncateFactor * vtxMatrices(vtx, 0) * dual_triangle_area(vtx, 0)/(radius*radius); + + double X = vtxCoords(vtx, 0)/radius; + double Y = vtxCoords(vtx, 1)/radius; + double Z = vtxCoords(vtx, 2)/radius; + if(isRotated){ + X = -vtxCoords(vtx, 2)/radius; + Z = vtxCoords(vtx, 0)/radius; + } - Kokkos::View m11_d("m11", vtxPerElm, nCells); - Kokkos::View m12_d("m12", vtxPerElm, nCells); - Kokkos::View m13_d("m13", vtxPerElm, nCells); - Kokkos::View m14_d("m14", vtxPerElm, nCells); - Kokkos::View m22_d("m22", vtxPerElm, nCells); - Kokkos::View m23_d("m23", vtxPerElm, nCells); - Kokkos::View m24_d("m23", vtxPerElm, nCells); - Kokkos::View m33_d("m33", vtxPerElm, nCells); - Kokkos::View m34_d("m34", vtxPerElm, nCells); - Kokkos::View m44_d("m34", vtxPerElm, nCells); - - auto sub_assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { - if(mask && (elm2Process(elm)==comm_rank)) { //if material point is 'active'/'enabled' - int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element - for(int i=0; iparallel_for(sub_assemble, "sub_assembly"); - pumipic::RecordTime("VtxSubAssemblyComputeCoeff", timer.seconds()); - - Kokkos::Timer timer2; - kkDbl2dViewHostU m11_h(m11, vtxPerElm, nCells); - kkDbl2dViewHostU m12_h(m12, vtxPerElm, nCells); - kkDbl2dViewHostU m13_h(m13, vtxPerElm, nCells); - kkDbl2dViewHostU m14_h(m14, vtxPerElm, nCells); - kkDbl2dViewHostU m22_h(m22, vtxPerElm, nCells); - kkDbl2dViewHostU m23_h(m23, vtxPerElm, nCells); - kkDbl2dViewHostU m24_h(m24, vtxPerElm, nCells); - kkDbl2dViewHostU m33_h(m33, vtxPerElm, nCells); - kkDbl2dViewHostU m34_h(m34, vtxPerElm, nCells); - kkDbl2dViewHostU m44_h(m44, vtxPerElm, nCells); - - Kokkos::deep_copy(m11_h, m11_d); - Kokkos::deep_copy(m12_h, m12_d); - Kokkos::deep_copy(m13_h, m13_d); - Kokkos::deep_copy(m14_h, m14_d); - Kokkos::deep_copy(m22_h, m22_d); - Kokkos::deep_copy(m23_h, m23_d); - Kokkos::deep_copy(m24_h, m24_d); - Kokkos::deep_copy(m33_h, m33_d); - Kokkos::deep_copy(m34_h, m34_d); - Kokkos::deep_copy(m44_h, m44_d); - pumipic::RecordTime("VtxSubAssemblyGetCoeff", timer2.seconds()); - -} -//Method 2 for coefficients Solve matrix -void MPMesh::solveMatrixAndRegularize(int nVertices, double* m11, double* m12, double* m13, double* m14, - double* m22, double* m23, double* m24, - double* m33, double* m34, - double* m44){ + auto diffTr = subM1(0, 0)-subM1(1, 1); + auto delG = sqrt(4.0 * pow(subM1(0, 1), 2) + pow(diffTr, 2)); + if(delG m11_h(m11, nVertices); - kkViewHostU m12_h(m12, nVertices); - kkViewHostU m13_h(m13, nVertices); - kkViewHostU m14_h(m14, nVertices); - kkViewHostU m22_h(m22, nVertices); - kkViewHostU m23_h(m23, nVertices); - kkViewHostU m24_h(m24, nVertices); - kkViewHostU m33_h(m33, nVertices); - kkViewHostU m34_h(m34, nVertices); - kkViewHostU m44_h(m44, nVertices); - - Kokkos::View m11_d("m11", nVertices); - Kokkos::View m12_d("m12", nVertices); - Kokkos::View m13_d("m13", nVertices); - Kokkos::View m14_d("m14", nVertices); - Kokkos::View m22_d("m22", nVertices); - Kokkos::View m23_d("m23", nVertices); - Kokkos::View m24_d("m24", nVertices); - Kokkos::View m33_d("m33", nVertices); - Kokkos::View m34_d("m34", nVertices); - Kokkos::View m44_d("m44", nVertices); - - Kokkos::deep_copy(m11_d, m11_h); - Kokkos::deep_copy(m12_d, m12_h); - Kokkos::deep_copy(m13_d, m13_h); - Kokkos::deep_copy(m14_d, m14_h); - Kokkos::deep_copy(m22_d, m22_h); - Kokkos::deep_copy(m23_d, m23_h); - Kokkos::deep_copy(m24_d, m24_h); - Kokkos::deep_copy(m33_d, m33_h); - Kokkos::deep_copy(m34_d, m34_h); - Kokkos::deep_copy(m44_d, m44_h); - pumipic::RecordTime("polyMPOsolveMatrixCoeffSet", timer.seconds()); - - Kokkos::Timer timer2; - auto dual_triangle_area=p_mesh->getMeshField(); - Kokkos::View VtxCoeffs("VtxCoeffs", nVertices); - double radius=p_mesh->getSphereRadius(); - Kokkos::parallel_for("fill", nVertices, KOKKOS_LAMBDA(const int vtx){ - Vec4d v0 = {m11_d(vtx), m12_d(vtx), m13_d(vtx), m14_d(vtx)}; - Vec4d v1 = {m12_d(vtx), m22_d(vtx), m23_d(vtx), m24_d(vtx)}; - Vec4d v2 = {m13_d(vtx), m23_d(vtx), m33_d(vtx), m34_d(vtx)}; - Vec4d v3 = {m14_d(vtx), m24_d(vtx), m34_d(vtx), m44_d(vtx)}; - //Matrix4d A = {v0,v1,v2,v3}; - Matrix4d A_regularized = {v0, v1, v2, v3}; - double regParam = sqrt(EPSILON)*(m11_d(vtx) + m22_d(vtx) + m33_d(vtx) + m44_d(vtx)); - A_regularized.addToDiag(regParam); - - double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0}; - CholeskySolve4d_UnitRHS(A_regularized, coeff); - - double mScale=sqrt(dual_triangle_area(vtx,0))/radius; - coeff[0]=coeff[0]*mScale*mScale; - coeff[1]=coeff[1]*mScale; - coeff[2]=coeff[2]*mScale; - coeff[3]=coeff[3]*mScale; - - for (int i=0; iprecomputedVtxCoeffs = VtxCoeffs; - pumipic::RecordTime("polyMPOsolveMatrixCoeffCompute", timer2.seconds()); - + this->precomputedVtxCoeffs_new = VtxCoeffs; } -//Method2 template -void MPMesh::subAssemblyVtx1(int vtxPerElm, int nCells, int comp, double* array) { - Kokkos::Timer timer; - - auto VtxCoeffs=this->precomputedVtxCoeffs; - - // MPI Information - MPI_Comm comm = p_MPs->getMPIComm(); - int comm_rank; - MPI_Comm_rank(comm, &comm_rank); - - //Mesh Information - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - int numVtx = p_mesh->getNumVertices(); - auto vtxCoords = p_mesh->getMeshField(); - auto elm2Process = p_mesh->getElm2Process(); - - // Material Points Information - constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice; - auto mpData = p_MPs->getData(); - auto weight = p_MPs->getData(); - auto mpPositions = p_MPs->getData(); - - double radius=p_mesh->getSphereRadius(); +void MPMesh::assemblyVtx1(){ + std::cout<<__FUNCTION__< array_d("reconstructedIceArea", vtxPerElm, nCells); - auto sub_assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { - if(mask && (elm2Process(elm)==comm_rank)) { - int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element - for(int i=0; igetMPIComm(); + MPI_Comm_rank(comm, &self); + MPI_Comm_size(comm, &numProcsTot); - auto factor = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] + - VtxCoeffs(vID,2)*CoordDiffs[2] + - VtxCoeffs(vID,3)*CoordDiffs[3]); - - auto val = factor*mpData(mp, comp); - Kokkos::atomic_add(&array_d(i, elm), val); - } - } - }; - p_MPs->parallel_for(sub_assemble, "sub_assembly"); - pumipic::RecordTime("polyMPOsubAssemblyFieldCompute", timer.seconds()); - - Kokkos::Timer timer2; - kkDbl2dViewHostU arrayHost(array, vtxPerElm, nCells); - Kokkos::deep_copy(arrayHost, array_d); - pumipic::RecordTime("PolyMPOsubAssemblyFieldGet", timer2.seconds()); -} + auto VtxCoeffs_new=this->precomputedVtxCoeffs_new; -//Method 1 -template -void MPMesh::assemblyVtx1() { - Kokkos::Timer timer; - //If no reconstruction till now calculate the coeffs - if (!isPreComputed) { - computeMatricesAndSolve(); - isPreComputed=true; - } - - auto VtxCoeffs=this->precomputedVtxCoeffs; //Mesh Information - auto elm2VtxConn = p_mesh->getElm2VtxConn(); + auto elm2VtxConn = p_mesh->getElm2VtxConn(); int numVtx = p_mesh->getNumVertices(); auto vtxCoords = p_mesh->getMeshField(); + int numVertices = p_mesh->getNumVertices(); //Mesh Field constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice; @@ -463,9 +315,6 @@ void MPMesh::assemblyVtx1() { if(p_mesh->getGeomType() == geom_spherical_surf) radius=p_mesh->getSphereRadius(); - //Reconstructed values - Kokkos::View reconVals("meshField", p_mesh->getNumVertices(), numEntries); - //Reconstruct auto reconstruct = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { if(mask) { //if material point is 'active'/'enabled' @@ -473,29 +322,350 @@ void MPMesh::assemblyVtx1() { for(int i=0; iparallel_for(reconstruct, "reconstruct"); + Kokkos::fence(); + pumipic::RecordTime("Assemble Field per process" + std::to_string(self), timer.seconds()); + + timer.reset(); + if(numProcsTot>1) + communicate_and_take_halo_contributions(meshField, numVertices, numEntries, 0, 0); + pumipic::RecordTime("Communicate Field Values" + std::to_string(self), timer.seconds()); +} - //Assign as a field - Kokkos::parallel_for("assigning", numVtx, KOKKOS_LAMBDA(const int vtx){ - for(int k=0; kgetMPIComm(); + MPI_Comm_rank(comm, &self); + MPI_Comm_size(comm, &numProcsTot); + + //The routine should work for elements too, although currently the communication + //is done for vertices. For elements, the follwoing three variables should correspond + //to elements. + auto entOwners = p_mesh->getVtx2Process(); + auto ent2global = p_mesh->getVtxGlobal(); + int numEntities = p_mesh->getNumVertices(); + + //Loop over elements and find no of owners and halos + Kokkos::View owner_count("owner_count"); + Kokkos::View halo_count("halo_count"); + Kokkos::deep_copy(owner_count, 0); + Kokkos::deep_copy(halo_count, 0); + Kokkos::parallel_for("countOwnerHalo", numEntities, KOKKOS_LAMBDA(const int elm){ + if (entOwners(elm)==self) + Kokkos::atomic_add(&owner_count(), 1); + else + Kokkos::atomic_add(&halo_count(), 1); }); - pumipic::RecordTime("PolyMPO_Reconstruct_Vtx1", timer.seconds()); + + Kokkos::deep_copy(numOwnersTot, owner_count); + Kokkos::deep_copy(numHalosTot, halo_count); + assert(numHalosTot+numOwnersTot == numEntities); + int num_ints_per_copy = 2; + + //#Halo Cells/proc which are owners on other process + numOwnersOnOtherProcs.resize(numProcsTot); + + //#OwnerCells/proc which are halos on other proces + numHalosOnOtherProcs.resize(numProcsTot); + + //For every halo cell find the owning process and the local Id in that process + haloOwnerProcs.reserve(numHalosTot); + haloOwnerLocalIDs.resize(numProcsTot); + + //Copy owning processes and globalIds to CPU + auto entOwners_host = Kokkos::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace::memory_space(), + entOwners); + auto ent2global_host = Kokkos::create_mirror_view_and_copy(Kokkos::DefaultHostExecutionSpace::memory_space(), + ent2global); + + //Do Map of Global To Local ID + //Check unordered vs ordered map; which faster? + std::map global2local; + //std::unordered_map global2local; + for (int iEnt = 0; iEnt < numEntities; iEnt++) { + int globalID = ent2global_host(iEnt); + global2local[globalID] = iEnt; + } + + //Loop over all halo Entities and find the owning process + for (auto iEnt=numOwnersTot; iEnt> sendBufs(numProcsTot); + for (int proc = 0; proc < numProcsTot; proc++) + sendBufs[proc].reserve(num_ints_per_copy*numOwnersOnOtherProcs[proc]); + + for (int iEnt=numOwnersTot; iEnt requests; + requests.reserve(2*numProcsTot); + + //Receive Calls + std::vector> recvBufs(numProcsTot); + for (int proc = 0; proc < numProcsTot; proc++) { + if (numHalosOnOtherProcs[proc] > 0) { + recvBufs[proc].resize(num_ints_per_copy*numHalosOnOtherProcs[proc]); + MPI_Request req; + MPI_Irecv(recvBufs[proc].data(), num_ints_per_copy*numHalosOnOtherProcs[proc], MPI_INT, proc, MPI_ANY_TAG, comm, &req); + requests.push_back(req); + } + } + //Send Calls + for (int proc=0; proc 0) { + ownerOwnerLocalIDs[proc].resize(numHalosOnOtherProcs[proc]); + ownerHaloLocalIDs[proc].resize(numHalosOnOtherProcs[proc]); + for (int i = 0; i < numHalosOnOtherProcs[proc]; i++) { + int globalID = recvBufs[proc][i*num_ints_per_copy]; + ownerOwnerLocalIDs[proc][i] = global2local[globalID]; + ownerHaloLocalIDs[proc][i] = recvBufs[proc][i*num_ints_per_copy+1]; + } + } + } + + // On the halo side, need to receive the localIds of owning Process + for (int proc = 0; proc < numProcsTot; proc++) { + if (numOwnersOnOtherProcs[proc] > 0) { // these are cells whose owners are in other processes + haloOwnerLocalIDs[proc].resize(numOwnersOnOtherProcs[proc]); + MPI_Request req; + MPI_Irecv(haloOwnerLocalIDs[proc].data(), haloOwnerLocalIDs[proc].size(), MPI_INT, proc, MPI_ANY_TAG, comm, &req); + requests.push_back(req); + } + } + + //Sends back localID of the owned cells so that HaloToOwner can be done for halo processes + for (int proc = 0; proc < numProcsTot; proc++) { + if (numHalosOnOtherProcs[proc]>0) { + MPI_Request req; + MPI_Isend(ownerOwnerLocalIDs[proc].data(), ownerOwnerLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req); + requests.push_back(req); + } + } + + MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); + + pumipic::RecordTime("Start Communication" + std::to_string(self), timer.seconds()); +} + +void MPMesh::communicate_and_take_halo_contributions(const Kokkos::View& meshField, int nEntities, int numEntries, int mode, int op){ + int self; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); + + Kokkos::Timer timer; + auto reconVals_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), meshField); + Kokkos::fence(); + std::vector> fieldData(nEntities, std::vector(numEntries, 0.0)); + for (int i = 0; i < nEntities; ++i) { + for (int j = 0; j < numEntries; ++j) { + fieldData[i][j] = reconVals_host(i, j); + } + } + pumipic::RecordTime("Communication-GPU to CPU-E-" + std::to_string(numEntries) + "-" + std::to_string(self), timer.seconds()); + + timer.reset(); + std::vector> recvIDVec; + std::vector> recvDataVec; + communicateFields(fieldData, nEntities, numEntries, mode, recvIDVec, recvDataVec); + pumipic::RecordTime("Communication-InterProcess-E-" + std::to_string(numEntries) + "-" + std::to_string(self), timer.seconds()); + + timer.reset(); + int numProcsTot = recvIDVec.size(); + //Flatten IDs + int totalSize = 0; + std::vector offsets(numProcsTot, 0); + for(int i=0; i flatIDVec(totalSize, 0); + for(int i=0; i recvIDGPU("recvIDGPU", totalSize); + auto hostView = Kokkos::View("recvIDCPU", totalSize); + std::copy(flatIDVec.begin(), flatIDVec.end(), hostView.data()); + Kokkos::deep_copy(recvIDGPU, hostView); + + //Flatten Data + int totalSize_data=0; + std::vector offsets_data(numProcsTot, 0); + for(int i=0; i flatDataVec(totalSize_data, 0); + for(int i=0; i recvDataGPU("recvDataGPU", totalSize_data); + auto hostView_data= Kokkos::View("recvDataCPU", totalSize_data); + std::copy(flatDataVec.begin(), flatDataVec.end(), hostView_data.data()); + Kokkos::deep_copy(recvDataGPU, hostView_data); + Kokkos::fence(); + //Assertions + assert(totalSize_data == totalSize*numEntries); + for (int i=0; i>& fieldData, const int numEntities, const int numEntries, int mode, + std::vector>& recvIDVec, std::vector>& recvDataVec){ + int self, numProcsTot; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); + MPI_Comm_size(comm, &numProcsTot); + + assert(numEntities == numOwnersTot + numHalosTot); + + std::vector> sendDataVec(numProcsTot); + + recvIDVec.resize(numProcsTot); + recvDataVec.resize(numProcsTot); + + for(int i = 0; i < numProcsTot; i++){ + if(i==self) continue; + + int numToSend = 0, numToRecv = 0; + if(mode == 0) { + //gather (halos send to owners) + numToSend = numOwnersOnOtherProcs[i]; + numToRecv = numHalosOnOtherProcs[i]; + } + else{ + //scatter (owners send to halos) + numToSend = numHalosOnOtherProcs[i]; + numToRecv = numOwnersOnOtherProcs[i]; + } + + if(numToSend > 0){ + sendDataVec[i].reserve(numToSend*numEntries); + } + if(numToRecv > 0){ + recvDataVec[i].resize(numToRecv*numEntries); + recvIDVec[i].resize(numToRecv); + } + } + + if(mode == 0){ + // Halos sends to owners + for (int iEnt = 0; iEnt < numHalosTot; iEnt++){ + auto ownerProc = haloOwnerProcs[iEnt]; + for (int iDouble = 0; iDouble < numEntries; iDouble++) + sendDataVec[ownerProc].push_back(fieldData[numOwnersTot+iEnt][iDouble]); + } + } + else if(mode == 1){ + // Owner sends to halos + for (size_t iProc=0; iProc requests; + requests.reserve(4*numProcsTot); + for(int proc = 0; proc < numProcsTot; proc++){ + if(proc == self) continue; + if(mode == 0 && numHalosOnOtherProcs[proc]){ + assert(recvIDVec[proc].size() == (size_t)numHalosOnOtherProcs[proc]); + assert(recvDataVec[proc].size() == recvIDVec[proc].size() * (size_t)numEntries); + MPI_Request req3, req4; + MPI_Irecv(recvIDVec[proc].data(), recvIDVec[proc].size(), MPI_INT, proc, 1, comm, &req3); + MPI_Irecv(recvDataVec[proc].data(), recvDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req4); + requests.push_back(req3); + requests.push_back(req4); + } + if(mode == 0 && numOwnersOnOtherProcs[proc]) { + assert(haloOwnerLocalIDs[proc].size() == (size_t)numOwnersOnOtherProcs[proc]); + assert(sendDataVec[proc].size() == haloOwnerLocalIDs[proc].size() * (size_t)numEntries); + MPI_Request req1, req2; + MPI_Isend(haloOwnerLocalIDs[proc].data(), haloOwnerLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req1); + MPI_Isend(sendDataVec[proc].data(), sendDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req2); + requests.push_back(req1); + requests.push_back(req2); + } + + if(mode == 1 && numOwnersOnOtherProcs[proc]){ + MPI_Request req3, req4; + MPI_Irecv(recvIDVec[proc].data(), recvIDVec[proc].size(), MPI_INT, proc, 1, comm, &req3); + MPI_Irecv(recvDataVec[proc].data(), recvDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req4); + requests.push_back(req3); + requests.push_back(req4); + } + if(mode == 1 && numHalosOnOtherProcs[proc]) { + MPI_Request req1, req2; + MPI_Isend(ownerHaloLocalIDs[proc].data(), ownerHaloLocalIDs[proc].size(), MPI_INT, proc, 1, comm, &req1); + MPI_Isend(sendDataVec[proc].data(), sendDataVec[proc].size(), MPI_DOUBLE, proc, 2, comm, &req2); + requests.push_back(req1); + requests.push_back(req2); + } + } + + MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } template @@ -519,98 +689,97 @@ void MPMesh::assembly(int order, MeshFieldType type, bool basisWeightFlag, bool // (HDT) weighted assembly of scalar field template DoubleView MPMesh::wtScaAssembly(){ - auto vtxCoords = p_mesh->getMeshField(); - int numVtxs = p_mesh->getNumVertices(); // total number of vertices of the mesh - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - auto mpPositions = p_MPs->getData(); - - DoubleView vField("wtScaField", numVtxs); // Kokkos array of double type, size = numVtxs - - auto mpData = p_MPs->getData(); - - auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { - if (mask) { - /* get the coordinates of all the vertices of elm */ - int nElmVtxs = elm2VtxConn(elm,0); // number of vertices bounding the element - Vec2d eVtxCoords[maxVtxsPerElm + 1]; - for (int i = 1; i <= nElmVtxs; i++) { - // elm2VtxConn(elm,i) is the vertex ID (1-based index) of vertex #i of elm - eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); - eVtxCoords[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); - } - // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) - eVtxCoords[nElmVtxs][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); - eVtxCoords[nElmVtxs][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); - - /* compute the values of basis functions at mp position */ - double basisByArea[maxElmsPerVtx]; - Vec2d mpCoord(mpPositions(mp,0), mpPositions(mp,1)); - getBasisByAreaGblForm(mpCoord, nElmVtxs, eVtxCoords, basisByArea); - - /* get the mp's property that is assembled to vertices */ - double assemVal = mpData(mp, 0); // ??? for scalar mp data, is index 0 always? - - /* accumulate the mp's property to vertices */ - for (int i = 0; i < nElmVtxs; i++) { - int vID = elm2VtxConn(elm,i+1)-1; - Kokkos::atomic_add(&vField(vID), assemVal * basisByArea[i]); - } + auto vtxCoords = p_mesh->getMeshField(); + int numVtxs = p_mesh->getNumVertices(); // total number of vertices of the mesh + auto elm2VtxConn = p_mesh->getElm2VtxConn(); + auto mpPositions = p_MPs->getData(); + + DoubleView vField("wtScaField", numVtxs); // Kokkos array of double type, size = numVtxs + + auto mpData = p_MPs->getData(); + + auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { + if (mask) { + /* get the coordinates of all the vertices of elm */ + int nElmVtxs = elm2VtxConn(elm,0); // number of vertices bounding the element + Vec2d eVtxCoords[maxVtxsPerElm + 1]; + for (int i = 1; i <= nElmVtxs; i++) { + // elm2VtxConn(elm,i) is the vertex ID (1-based index) of vertex #i of elm + eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); + eVtxCoords[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); + } + // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) + eVtxCoords[nElmVtxs][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); + eVtxCoords[nElmVtxs][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); + + /* compute the values of basis functions at mp position */ + double basisByArea[maxElmsPerVtx]; + Vec2d mpCoord(mpPositions(mp,0), mpPositions(mp,1)); + getBasisByAreaGblForm(mpCoord, nElmVtxs, eVtxCoords, basisByArea); + + /* get the mp's property that is assembled to vertices */ + double assemVal = mpData(mp, 0); // ??? for scalar mp data, is index 0 always? + + /* accumulate the mp's property to vertices */ + for (int i = 0; i < nElmVtxs; i++) { + int vID = elm2VtxConn(elm,i+1)-1; + Kokkos::atomic_add(&vField(vID), assemVal * basisByArea[i]); } - }; - p_MPs->parallel_for(assemble, "weightedScalarAssembly"); - return vField; + } + }; + p_MPs->parallel_for(assemble, "weightedScalarAssembly"); + return vField; } // wtScaAssembly - // (HDT) weighted assembly of vector2 field (not weighted by mass/volume) template Vec2dView MPMesh::wtVec2Assembly(){ - auto vtxCoords = p_mesh->getMeshField(); - int numVtxs = p_mesh->getNumVertices(); // total number of vertices of the mesh - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - auto mpPositions = p_MPs->getData(); - - Vec2dView vField("wtVec2Field", numVtxs); // Kokkos array of Vec2d type, size = numVtxs - - auto mpData = p_MPs->getData(); - - auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { - if (mask) { - /* collect the coordinates of all the vertices of elm */ - int nElmVtxs = elm2VtxConn(elm,0); // number of vertices bounding the element - Vec2d eVtxCoords[maxVtxsPerElm + 1]; - for (int i = 1; i <= nElmVtxs; i++) { - // elm2VtxConn(elm,i) is the vertex ID (1-based index) of vertex #i of elm - eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); - eVtxCoords[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); - } - // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) - eVtxCoords[nElmVtxs][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); - eVtxCoords[nElmVtxs][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); - - /* compute the values of basis functions at mp position */ - double basisByArea[maxElmsPerVtx]; - Vec2d mpCoord(mpPositions(mp,0), mpPositions(mp,1)); - getBasisByAreaGblForm(mpCoord, nElmVtxs, eVtxCoords, basisByArea); - - /* get the mp's volume */ - double mpVolume = 1.0; // TODO: change to mp's volume here - - /* get the mp's property to be assembled */ - Vec2d assemVal; - assemVal[0] = mpData(mp, 0) * mpVolume; - assemVal[1] = mpData(mp, 1) * mpVolume; - - /* accumulate the mp's constructed quantities to the cell vertices */ - for (int i = 0; i < nElmVtxs; i++) { - int vID = elm2VtxConn(elm,i+1)-1; - Kokkos::atomic_add(&(vField(vID)[0]), assemVal[0] * basisByArea[i]); - Kokkos::atomic_add(&(vField(vID)[1]), assemVal[1] * basisByArea[i]); - } + auto vtxCoords = p_mesh->getMeshField(); + int numVtxs = p_mesh->getNumVertices(); // total number of vertices of the mesh + auto elm2VtxConn = p_mesh->getElm2VtxConn(); + auto mpPositions = p_MPs->getData(); + + Vec2dView vField("wtVec2Field", numVtxs); // Kokkos array of Vec2d type, size = numVtxs + + auto mpData = p_MPs->getData(); + + auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { + if (mask) { + /* collect the coordinates of all the vertices of elm */ + int nElmVtxs = elm2VtxConn(elm,0); // number of vertices bounding the element + Vec2d eVtxCoords[maxVtxsPerElm + 1]; + for (int i = 1; i <= nElmVtxs; i++) { + // elm2VtxConn(elm,i) is the vertex ID (1-based index) of vertex #i of elm + eVtxCoords[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); + eVtxCoords[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); } - }; - p_MPs->parallel_for(assemble, "weightedVec2dAssembly"); - return vField; + // last component of eVtxCoords stores the firs vertex (to avoid if-condition in the Wachspress computation) + eVtxCoords[nElmVtxs][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); + eVtxCoords[nElmVtxs][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); + + /* compute the values of basis functions at mp position */ + double basisByArea[maxElmsPerVtx]; + Vec2d mpCoord(mpPositions(mp,0), mpPositions(mp,1)); + getBasisByAreaGblForm(mpCoord, nElmVtxs, eVtxCoords, basisByArea); + + /* get the mp's volume */ + double mpVolume = 1.0; // TODO: change to mp's volume here + + /* get the mp's property to be assembled */ + Vec2d assemVal; + assemVal[0] = mpData(mp, 0) * mpVolume; + assemVal[1] = mpData(mp, 1) * mpVolume; + + /* accumulate the mp's constructed quantities to the cell vertices */ + for (int i = 0; i < nElmVtxs; i++) { + int vID = elm2VtxConn(elm,i+1)-1; + Kokkos::atomic_add(&(vField(vID)[0]), assemVal[0] * basisByArea[i]); + Kokkos::atomic_add(&(vField(vID)[1]), assemVal[1] * basisByArea[i]); + } + } + }; + p_MPs->parallel_for(assemble, "weightedVec2dAssembly"); + return vField; } // wtVec2Assembly template diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index 2369d57..0f15b92 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -13,6 +13,7 @@ namespace{ } } +//initialize and finalize void polympo_initialize_f() { int isMPIInit; MPI_Initialized(&isMPIInit); @@ -24,6 +25,7 @@ void polympo_finalize_f() { Kokkos::finalize(); } +//create/delete MpMesh object MPMesh_ptr polympo_createMPMesh_f(const int testMeshOption, const int testMPOption) { polyMPO::Mesh* p_mesh; if(testMeshOption){ @@ -52,6 +54,7 @@ void polympo_deleteMPMesh_f(MPMesh_ptr p_mpmesh) { delete (polyMPO::MPMesh*)p_mpmesh; } +//set MPI communicator void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm){ checkMPMeshValid(p_mpmesh); MPI_Comm comm = MPI_Comm_f2c(fcomm); @@ -59,12 +62,18 @@ void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm){ p_MPs->setMPIComm(comm); } +void polympo_startCommunication_f(MPMesh_ptr p_mpmesh){ + auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); + mpmesh->startCommunication(); //Temporary last API in set mesh needs to be separate API +} + +//MP info void polympo_createMPs_f(MPMesh_ptr p_mpmesh, - const int numElms, - const int numMPs, // total number of MPs which is >= number of active MPs - int* mpsPerElm, - const int* mp2Elm, - const int* isMPActive) { + const int numElms, + const int numMPs, // total number of MPs which is >= number of active MPs + int* mpsPerElm, + const int* mp2Elm, + const int* isMPActive){ checkMPMeshValid(p_mpmesh); //the mesh must be fixed/set before adding MPs auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; @@ -87,7 +96,7 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh, MPI_Allreduce(&numActiveMPs, &globalNumActiveMPs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(&minElmID, &globalMinElmID, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); PMT_ALWAYS_ASSERT(globalNumActiveMPs>0); - + //Loop over all mesh elements 0,1,... and find the first element that has an associated MP int firstElmWithMPs=INT_MAX; for (int i=0; ip_MPs; - p_MPs->setElmIDoffset(offset); - + p_MPs->setElmIDoffset(offset); } void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, - const int numMPs, // Total MPs which is GREATER than or equal to number of active MPs - const int* allMP2Elm, - const int* addedMPMask) { + const int numMPs, // Total MPs which is GREATER than or equal to number of active MPs + const int* allMP2Elm, + const int* addedMPMask){ checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); @@ -196,13 +204,13 @@ void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, } void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh, - const int sizeMP2elm, - const int* elem_ids, - const int nMPs_delete, - const int nMPs_add, - int* recvMPs_elm, - int* recvMPs_ids) { - + const int sizeMP2elm, + const int* elem_ids, + const int nMPs_delete, + const int nMPs_add, + int* recvMPs_elm, + int* recvMPs_ids) { + Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -216,7 +224,7 @@ void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh, auto elem_ids_d = create_mirror_view_and_copy(elem_ids, sizeMP2elm); auto recvMPs_elm_d = create_mirror_view_and_copy(recvMPs_elm, nMPs_add); auto recvMPs_ids_d = create_mirror_view_and_copy(recvMPs_ids, nMPs_add); - + Kokkos::View mp2Elm("mp2Elm", p_MPs->getCapacity()); Kokkos::View numDeletedMPs_d("numDeletedMPs", 1); auto mpAppID = p_MPs->getData(); @@ -235,6 +243,11 @@ void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh, pumipic::RecordTime("polympo_startRebuildMPs2_f", timer.seconds()); } +int polympo_getMPCount_f(MPMesh_ptr p_mpmesh) { + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + return p_MPs->getCount(); +} + void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); @@ -250,9 +263,7 @@ void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appI p_MPs->setAppIDFunc(getNextAppID); } -void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, - const int numMPs, - int* elmIDs){ +void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -275,10 +286,7 @@ void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, pumipic::RecordTime("PolyMPO_getMPTgtElmID", timer.seconds()); } -void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, - const int numMPs, - int* elmIDs){ - +void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs){ checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); @@ -296,21 +304,17 @@ void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, } }; p_MPs->parallel_for(getElmId, "get mpCurElmID"); - Kokkos::deep_copy( arrayHost, mpCurElmIDCopy); - + Kokkos::deep_copy( arrayHost, mpCurElmIDCopy); } void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag){ - //chech validity checkMPMeshValid(p_mpmesh); - ((polyMPO::MPMesh*)p_mpmesh)->p_MPs->setRotatedFlag(isRotateFlag>0); - + ((polyMPO::MPMesh*)p_mpmesh)->p_mesh->setRotatedFlag(isRotateFlag>0); } -void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, - const int nComps, - const int numMPs, - const double* mpPositionsIn){ +//MP Slices +//Positions +void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -340,10 +344,7 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, pumipic::RecordTime("PolyMPO_setMPPositions", timer.seconds()); } -void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh, - const int nComps, - const int numMPs, - double* mpPositionsHost){ +void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpPositionsHost){ checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries); @@ -365,11 +366,7 @@ void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh, Kokkos::deep_copy(arrayHost, mpPositionsCopy); } - -void polympo_setMPTgtPositions_f(MPMesh_ptr p_mpmesh, - const int nComps, - const int numMPs, - const double* mpPositionsIn){ +void polympo_setMPTgtPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -398,10 +395,7 @@ void polympo_setMPTgtPositions_f(MPMesh_ptr p_mpmesh, pumipic::RecordTime("PolyMPO_setMPTgtPositions", timer.seconds()); } -void polympo_getMPTgtPositions_f(MPMesh_ptr p_mpmesh, - const int nComps, - const int numMPs, - double* mpPositionsHost){ +void polympo_getMPTgtPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpPositionsHost){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -425,11 +419,8 @@ void polympo_getMPTgtPositions_f(MPMesh_ptr p_mpmesh, pumipic::RecordTime("PolyMPO_getMPTgtPositions", timer.seconds()); } - -void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, - const int nComps, - const int numMPs, - const double* mpRotLatLonIn){ +//LatLon +void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn){ Kokkos::Timer timer; static int callCount = 0; PMT_ALWAYS_ASSERT(callCount == 0); @@ -455,10 +446,7 @@ void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, pumipic::RecordTime("PolyMPO_setMPRotLatLon", timer.seconds()); } -void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, - const int nComps, - const int numMPs, - double* mpRotLatLonHost){ +void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost){ checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); @@ -479,11 +467,7 @@ void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, Kokkos::deep_copy(arrayHost, mpRotLatLonCopy); } - -void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, - const int nComps, - const int numMPs, - const double* mpRotLatLonIn){ +void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -506,10 +490,7 @@ void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, pumipic::RecordTime("PolyMPO_setMPTgtRotLatLon", timer.seconds()); } -void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, - const int nComps, - const int numMPs, - double* mpRotLatLonHost){ +void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost){ Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; @@ -532,11 +513,15 @@ void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, pumipic::RecordTime("PolyMPO_getMPTgtRotLatLon", timer.seconds()); } - +//MP Fields void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn) { Kokkos::Timer timer; checkMPMeshValid(p_mpmesh); auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + int self; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); + PMT_ALWAYS_ASSERT(nComps == 1); //TODO mp_sclr_t PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); @@ -552,7 +537,7 @@ void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs } }; p_MPs->parallel_for(setMPMass, "setMPMass"); - pumipic::RecordTime("PolyMPO_setMPMass", timer.seconds()); + pumipic::RecordTime("PolyMPO_setMPMass" + std::to_string(self), timer.seconds()); } void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpMassHost) { @@ -623,31 +608,12 @@ void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, pumipic::RecordTime("PolyMPO_getMPVel", timer.seconds()); } -//TODO: implement these -void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStrainRateIn){ +void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); - auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; - PMT_ALWAYS_ASSERT(nComps == 6); //TODO: mp_sym_mat3d_t - PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount()); - PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID()); - - auto mpStrainRate = p_MPs->getData(); - auto mpAppID = p_MPs->getData(); - kkViewHostU mpStrainRateIn_h(mpStrainRateIn,nComps,numMPs); - Kokkos::View mpStrainRateIn_d("mpStrainRateDevice",nComps,numMPs); - Kokkos::deep_copy(mpStrainRateIn_d, mpStrainRateIn_h); - auto setMPStrainRate = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ - if(mask){ - mpStrainRate(mp,0) = mpStrainRateIn_d(0, mpAppID(mp)); - mpStrainRate(mp,1) = mpStrainRateIn_d(1, mpAppID(mp)); - mpStrainRate(mp,2) = mpStrainRateIn_d(2, mpAppID(mp)); - mpStrainRate(mp,3) = mpStrainRateIn_d(3, mpAppID(mp)); - mpStrainRate(mp,4) = mpStrainRateIn_d(4, mpAppID(mp)); - mpStrainRate(mp,5) = mpStrainRateIn_d(5, mpAppID(mp)); - } - }; - p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate"); + auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh); + mpMesh->calculateStrain(); } + void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStrainRateHost){ checkMPMeshValid(p_mpmesh); std::cerr << "Error: This routine is not implemented yet\n"; @@ -657,6 +623,7 @@ void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int (void)numMPs; (void)mpStrainRateHost; } + void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn){ checkMPMeshValid(p_mpmesh); std::cerr << "Error: This routine is not implemented yet\n"; @@ -666,6 +633,7 @@ void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numM (void)numMPs; (void)mpStressIn; } + void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost){ checkMPMeshValid(p_mpmesh); std::cerr << "Error: This routine is not implemented yet\n"; @@ -817,6 +785,64 @@ void polympo_setMeshElm2ElmConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const }); } +//Owning Process and Global IDs +void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array){ + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); + kkViewHostU arrayHost(array,nCells); + + //check the size + PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); + + Kokkos::View owningProc("owningProc",nCells); + Kokkos::deep_copy(owningProc, arrayHost); + p_mesh->setOwningProc(owningProc); +} + +void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + kkViewHostU arrayHost(array,nVertices); + + //check the size + PMT_ALWAYS_ASSERT(nVertices == p_mesh->getNumVertices()); + + Kokkos::View owningProcVertex("owningProcVerterx", nVertices); + Kokkos::deep_copy(owningProcVertex, arrayHost); + p_mesh->setOwningProcVertex(owningProcVertex); +} + +void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array){ + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + Kokkos::View arrayHost("arrayHost", nCells); + for (int i = 0; i < nCells; i++) { + arrayHost(i) = array[i] - 1; // TODO right now elmID offset is set after MPs initialized + } + //check the size + PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); + + Kokkos::View elmGlobal("elmGlobal",nCells); + Kokkos::deep_copy(elmGlobal, arrayHost); + p_mesh->setElmGlobal(elmGlobal); +} + +void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + PMT_ALWAYS_ASSERT(nVertices==p_mesh->getNumVertices()); + Kokkos::View arrayHost("arrayHost", nVertices); + for (int i = 0; i < nVertices; i++) { + arrayHost(i) = array[i] - 1; // TODO right now elmID offset is set after MPs initialized + } + + Kokkos::View vtxGlobal("vtxGlobal",nVertices); + Kokkos::deep_copy(vtxGlobal, arrayHost); + p_mesh->setVtxGlobal(vtxGlobal); +} + +//Mesh Fields int polympo_getMeshFVtxType_f() { return polyMPO::MeshFType_VtxBased; } @@ -901,75 +927,6 @@ void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double } } -void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const double* xArray, const double* yArray, const double* zArray){ - //chech validity - checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - - //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); - - //copy the host array to the device - auto elmCenter = p_mesh->getMeshField(); - auto h_elmCenter = Kokkos::create_mirror_view(elmCenter); - for(int i=0; ip_mesh; - - //check the size - PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); - - //copy the device to host - auto elmCenter = p_mesh->getMeshField(); - auto h_elmCenter = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), elmCenter); - for(int i=0; ip_mesh; - - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - //copy the host array to the device - auto dualArea = p_mesh->getMeshField(); - auto h_dualArea = Kokkos::create_mirror_view(dualArea); - for(int i=0; ip_mesh; - - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); - //copy the device to host - auto dualArea = p_mesh->getMeshField(); - auto h_dualArea = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dualArea); - for(int i=0; ip_mesh; - + auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; + int self; + MPI_Comm comm = p_MPs->getMPIComm(); + MPI_Comm_rank(comm, &self); + //check the size PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); @@ -1039,7 +1000,7 @@ void polympo_getMeshVtxMass_f(MPMesh_ptr p_mpmesh, const int nVertices, double* for(int i=0; ip_mesh; + kkDbl2dViewHostU arrayHost(array,nComps,nVertices); + Kokkos::View array_d("meshVelIncrDevice",nComps,nVertices); + + auto vtxField = p_mesh->getMeshField(); + + //check the size + PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); + PMT_ALWAYS_ASSERT(static_cast(nVertices*vec2d_nEntries)==vtxField.size()); + + //copy the device array to the host + Kokkos::parallel_for("get mesh dispIncr", nVertices, KOKKOS_LAMBDA(const int iVtx){ + array_d(0,iVtx) = vtxField(iVtx,0); + array_d(1,iVtx) = vtxField(iVtx,1); + }); + Kokkos::deep_copy(arrayHost, array_d); +} + void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) { Kokkos::Timer timer; //check mpMesh is valid @@ -1126,15 +1109,14 @@ void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c pumipic::RecordTime("PolyMPO_setMeshVtxOnSurfDispIncr", timer.seconds()); } - -void polympo_getMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array) { +void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array) { //check mpMesh is valid checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; kkDbl2dViewHostU arrayHost(array,nComps,nVertices); - Kokkos::View array_d("meshVelIncrDevice",nComps,nVertices); + Kokkos::View array_d("meshDispIncrDevice",nComps,nVertices); - auto vtxField = p_mesh->getMeshField(); + auto vtxField = p_mesh->getMeshField(); //check the size PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); @@ -1149,32 +1131,95 @@ void polympo_getMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c Kokkos::deep_copy(arrayHost, array_d); } -void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array) { - //check mpMesh is valid +void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const double* xArray, const double* yArray, const double* zArray){ + //chech validity checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - kkDbl2dViewHostU arrayHost(array,nComps,nVertices); - Kokkos::View array_d("meshDispIncrDevice",nComps,nVertices); - auto vtxField = p_mesh->getMeshField(); + //check the size + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); + + //copy the host array to the device + auto elmCenter = p_mesh->getMeshField(); + auto h_elmCenter = Kokkos::create_mirror_view(elmCenter); + for(int i=0; ip_mesh; //check the size - PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries); - PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); - PMT_ALWAYS_ASSERT(static_cast(nVertices*vec2d_nEntries)==vtxField.size()); + PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); - //copy the device array to the host - Kokkos::parallel_for("get mesh dispIncr", nVertices, KOKKOS_LAMBDA(const int iVtx){ - array_d(0,iVtx) = vtxField(iVtx,0); - array_d(1,iVtx) = vtxField(iVtx,1); - }); - Kokkos::deep_copy(arrayHost, array_d); + //copy the device to host + auto elmCenter = p_mesh->getMeshField(); + auto h_elmCenter = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), elmCenter); + for(int i=0; ipush1P(); - return is_migrating; + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + //copy the host array to the device + auto dualArea = p_mesh->getMeshField(); + auto h_dualArea = Kokkos::create_mirror_view(dualArea); + for(int i=0; ip_mesh; + + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + //copy the device to host + auto dualArea = p_mesh->getMeshField(); + auto h_dualArea = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dualArea); + for(int i=0; ip_mesh; + bool isRotated = p_mesh->getRotatedFlag(); + p_mesh->setGnomonicProjection(isRotated); +} + +void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){ + //chech validity + checkMPMeshValid(p_mpmesh); + auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; + + PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices); + //copy the host array to the device + auto tanLatVertexRotatedOverRadius = p_mesh->getMeshField(); + auto h_tanLatVertexRotatedOverRadius = Kokkos::create_mirror_view(tanLatVertexRotatedOverRadius); + for(int i=0; ipush(); } void polympo_push_ahead_f(MPMesh_ptr p_mpmesh){ @@ -1182,6 +1227,12 @@ void polympo_push_ahead_f(MPMesh_ptr p_mpmesh){ ((polyMPO::MPMesh*)p_mpmesh)->push_ahead(); } +bool polympo_push1P_f(MPMesh_ptr p_mpmesh){ + checkMPMeshValid(p_mpmesh); + bool is_migrating=((polyMPO::MPMesh*)p_mpmesh)->push1P(); + return is_migrating; +} + void polympo_push_swap_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); ((polyMPO::MPMesh*)p_mpmesh)->push_swap(); @@ -1192,12 +1243,7 @@ void polympo_push_swap_pos_f(MPMesh_ptr p_mpmesh){ ((polyMPO::MPMesh*)p_mpmesh)->push_swap_pos(); } -void polympo_push_f(MPMesh_ptr p_mpmesh){ - checkMPMeshValid(p_mpmesh); - ((polyMPO::MPMesh*)p_mpmesh)->push(); -} - -//TODO skeleton of reconstruction functions +// Reconstruction of variables from MPs to mesh vertices void polympo_setReconstructionOfMass_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType){ checkMPMeshValid(p_mpmesh); auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); @@ -1236,104 +1282,32 @@ void polympo_setReconstructionOfStress_f(MPMesh_ptr p_mpmesh, const int order, c (void)meshEntType; } -//With MPI communication done via MPAS -void polympo_vtxSubAssemblyIceArea_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int comp, double* array){ - checkMPMeshValid(p_mpmesh); - auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - PMT_ALWAYS_ASSERT(vtxPerElm <= maxVtxsPerElm); - PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); - PMT_ALWAYS_ASSERT(comp == 0 || comp== 1); //either first or second component - mpmesh->subAssemblyVtx1(vtxPerElm, nCells, comp, array); -} - -void polympo_vtxSubAssemblyVelocity_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int comp, double* array){ +void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - PMT_ALWAYS_ASSERT(vtxPerElm <= maxVtxsPerElm); - PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); - PMT_ALWAYS_ASSERT(comp == 0 || comp== 1); //either first or second component - mpmesh->subAssemblyVtx1(vtxPerElm, nCells, comp, array); + mpmesh->reconstructSlices(); } -void polympo_subAssemblyCoeffs_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, double* m11, double* m12, double* m13, double* m14, - double* m22, double* m23, double* m24, - double* m33, double* m34, - double* m44){ +//Simpler/cleaner way of calling reconstruction +void polympo_reconstruct_coeff_with_MPI_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - PMT_ALWAYS_ASSERT(vtxPerElm <= maxVtxsPerElm); - PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); - mpmesh->subAssemblyCoeffs(vtxPerElm, nCells, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44); + mpmesh->reconstruct_coeff_full(); } -void polympo_regularize_and_solve_matrix_f(MPMesh_ptr p_mpmesh, int nVertices, double* m11, double* m12, double* m13, double* m14, - double* m22, double* m23, double* m24, - double* m33, double* m34, - double* m44){ +void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - PMT_ALWAYS_ASSERT(nVertices == p_mesh->getNumVertices()); auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); - mpmesh->solveMatrixAndRegularize(nVertices, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44); + mpmesh->assemblyVtx1(); } -void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh){ +void polympo_reconstruct_velocity_with_MPI_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh); - mpmesh->reconstructSlices(); -} - -void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array){ - checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - PMT_ALWAYS_ASSERT(p_mesh->meshEditable()); - kkViewHostU arrayHost(array,nCells); - - //check the size - PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); - - Kokkos::View owningProc("owningProc",nCells); - Kokkos::deep_copy(owningProc, arrayHost); - p_mesh->setOwningProc(owningProc); -} - -void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array){ - checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - Kokkos::View arrayHost("arrayHost", nCells); - for (int i = 0; i < nCells; i++) { - arrayHost(i) = array[i] - 1; // TODO right now elmID offset is set after MPs initialized - } - //check the size - PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements()); - - Kokkos::View elmGlobal("elmGlobal",nCells); - Kokkos::deep_copy(elmGlobal, arrayHost); - p_mesh->setElmGlobal(elmGlobal); -} - -void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){ - checkMPMeshValid(p_mpmesh); - auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; - PMT_ALWAYS_ASSERT(nVertices==p_mesh->getNumVertices()); - Kokkos::View arrayHost("arrayHost", nVertices); - for (int i = 0; i < nVertices; i++) { - arrayHost(i) = array[i] - 1; // TODO right now elmID offset is set after MPs initialized - } - - Kokkos::View vtxGlobal("vtxGlobal",nVertices); - Kokkos::deep_copy(vtxGlobal, arrayHost); - p_mesh->setVtxGlobal(vtxGlobal); -} - -int polympo_getMPCount_f(MPMesh_ptr p_mpmesh) { - auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs; - return p_MPs->getCount(); + mpmesh->assemblyVtx1(); } +//Timing void polympo_enableTiming_f(){ pumipic::EnableTiming(); } diff --git a/src/pmpo_c.h b/src/pmpo_c.h index d39c195..ecca7dd 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -17,18 +17,15 @@ void polympo_deleteMPMesh_f(MPMesh_ptr p_mpmesh); //set MPI communicator void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm); //TODO: add a function to get communicator +void polympo_startCommunication_f(MPMesh_ptr p_mpmesh); //MP info void polympo_createMPs_f(MPMesh_ptr p_mpmesh, const int numElms, const int numMPs, int* mpsPerElm, const int* mp2Elm, const int* isMPActive); void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, const int numMPs, const int* allTgtMpElmIn, const int* addedMPMask); void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh, const int size1, const int* arg1, const int size2, const int size3, int* arg2, int* arg3); - int polympo_getMPCount_f(MPMesh_ptr p_mpmesh); - void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh); - void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs); - void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs); void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag); @@ -44,12 +41,12 @@ void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int n void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost); void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn); void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost); - +//MP fields void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn); void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpMassHost); void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn); void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost); -void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStrainRateIn); +void polympo_setMPStrainRate_f(MPMesh_ptr p_mpmesh); void polympo_getMPStrainRate_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStrainRateHost); void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn); void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost); @@ -71,10 +68,10 @@ void polympo_setMeshNumEdgesPerElm_f(MPMesh_ptr p_mpmesh, const int nCells, cons void polympo_setMeshElm2VtxConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int nCells, const int* array); void polympo_setMeshElm2ElmConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int nCells, const int* array); void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array); +void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array); void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array); void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array); - //Mesh fields int polympo_getMeshFVtxType_f(); int polympo_getMeshFElmType_f(); @@ -94,11 +91,12 @@ void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const double* xArray, const double* yArray, const double* zArray); void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* xArray, double* yArray, double* zArray); -//Area Triangle void polympo_setMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* areaTriangle); void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, double* areaTriangle); - -// calculations +void polympo_setGnomonicProjection_f(MPMesh_ptr p_mpmesh); +void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array); + +// Advection calculations void polympo_push_f(MPMesh_ptr p_mpmesh); void polympo_push_ahead_f(MPMesh_ptr p_mpmesh); bool polympo_push1P_f(MPMesh_ptr p_mpmesh); @@ -111,19 +109,11 @@ void polympo_setReconstructionOfVel_f(MPMesh_ptr p_mpmesh, const int order, cons void polympo_setReconstructionOfStrainRate_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType); void polympo_setReconstructionOfStress_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType); void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh); +//Simpler/cleaner way of calling reconstruction +void polympo_reconstruct_coeff_with_MPI_f(MPMesh_ptr p_mpmesh); +void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh); +void polympo_reconstruct_velocity_with_MPI_f(MPMesh_ptr p_mpmesh); -//Reconstruction using MPAS -void polympo_vtxSubAssemblyIceArea_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int comp, double* array); -void polympo_vtxSubAssemblyVelocity_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int comp, double* array); -void polympo_subAssemblyCoeffs_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, double* m11, double* m12, double* m13, double* m14, - double* m22, double* m23, double* m24, - double* m33, double* m34, - double* m44); -void polympo_regularize_and_solve_matrix_f(MPMesh_ptr p_mpmesh, int nVetices, double* m11, double* m12, double* m13, double* m14, - double* m22, double* m23, double* m24, - double* m33, double* m34, - double* m44); - // Timing void polympo_enableTiming_f(); void polympo_summarizeTime_f(); diff --git a/src/pmpo_createTestMPMesh.cpp b/src/pmpo_createTestMPMesh.cpp index 3a2ae6c..0f314bb 100644 --- a/src/pmpo_createTestMPMesh.cpp +++ b/src/pmpo_createTestMPMesh.cpp @@ -217,7 +217,7 @@ MaterialPoints* initTestMPs(Mesh* mesh, int testMPOption){ } if(geomType == geom_spherical_surf){ auto p_MPs = new MaterialPoints(numElms,numMPs,positions,numMPsPerElement,MPToElement); - p_MPs->setRotatedFlag(false); + mesh->setRotatedFlag(false); auto mpRotLatLonField = p_MPs->getData(); auto setRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){ mpRotLatLonField(mp,0) = latLonPositions(mp,0); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 10d04f2..2347717 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -13,6 +13,7 @@ module polympo subroutine polympo_initialize() bind(C, NAME='polympo_initialize_f') use :: iso_c_binding end subroutine + !--------------------------------------------------------------------------- !> @brief finalize polympo, no polympo apis may be called after this !> @remark the user must not finalize MPI until after this call @@ -20,6 +21,7 @@ subroutine polympo_initialize() bind(C, NAME='polympo_initialize_f') subroutine polympo_finalize() bind(C, NAME='polympo_finalize_f') use :: iso_c_binding end subroutine + !--------------------------------------------------------------------------- !> @brief create MPMesh object !> @param testMeshOption >0 init a planar test Mesh @@ -33,6 +35,7 @@ function polympo_createMPMesh(setMeshOption, setMPOption) bind(C, NAME='polympo_ type(c_ptr) polympo_createMPMesh integer(c_int), value :: setMeshOption, setMPOption end function + !--------------------------------------------------------------------------- !> @brief delete MPMesh object !> @param mpmesh(in/out) MPMesh object @@ -41,6 +44,7 @@ subroutine polympo_deleteMPMesh(mpMesh) bind(C, NAME='polympo_deleteMPMesh_f') use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief set the MPI communicator used by polympo !> @param comm(in) MPI communicator @@ -51,6 +55,13 @@ subroutine polympo_setMPICommunicator(mpMesh, comm) & type(c_ptr), value :: mpMesh integer(c_int), value :: comm end subroutine + + subroutine polympo_startCommunication(mpMesh) & + bind(C, NAME='polympo_startCommunication_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + !--------------------------------------------------------------------------- !> @brief create the material points !> @brief the fields associated with the MPs are NOT initialized @@ -71,6 +82,7 @@ subroutine polympo_createMPs(mpMesh, numElms, numMPs, mpsPerElm, mp2Elm, isMPAct type(c_ptr), intent(in), value :: mp2Elm type(c_ptr), intent(in), value :: isMPActive end subroutine + !--------------------------------------------------------------------------- !> @brief move MPs to a new element, add new MPs, or delete MPs !> @brief the fields associated with the MPs are NOT initialized @@ -88,7 +100,6 @@ subroutine polympo_startRebuildMPs(mpMesh, numMPs, allMP2Elm, addedMPMask) & type(c_ptr), intent(in), value :: addedMPMask end subroutine - subroutine polympo_startRebuildMPs2(mpMesh, size1, arg1, size2, size3, arg2, arg3) & bind(C, NAME='polympo_startRebuildMPs2_f') use :: iso_c_binding @@ -100,7 +111,13 @@ subroutine polympo_startRebuildMPs2(mpMesh, size1, arg1, size2, size3, arg2, arg type(c_ptr), value :: arg2 type(c_ptr), value :: arg3 end subroutine - ! + + integer function polympo_getMPCount(mpMesh) bind(C, name="polympo_getMPCount_f") + use :: iso_c_binding + implicit none + type(c_ptr), value :: mpMesh + end function + !--------------------------------------------------------------------------- !> @brief called after startRebuild() !> @brief called after initializing MP fields @@ -111,6 +128,7 @@ subroutine polympo_finishRebuildMPs(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief Stores pointer to appID data structure and a function to retrieve them used in migration !> @param mpmesh(in/out) MPMesh object @@ -124,7 +142,7 @@ subroutine polympo_setAppIDFunc(mpMesh, getNext, appIDs) & type(c_funptr), value :: getNext type(c_ptr), value :: appIDs end subroutine - + !--------------------------------------------------------------------------- !> @brief get the current element ID MP array from a polympo array !> @param mpmesh(in/out) MPMesh object @@ -139,7 +157,6 @@ subroutine polympo_getMPCurElmID(mpMesh, numMPs, array) & type(c_ptr), value :: array end subroutine - subroutine polympo_getMPTgtElmID(mpMesh, numMPs, array) & bind(C, NAME='polympo_getMPTgtElmID_f') use :: iso_c_binding @@ -147,7 +164,7 @@ subroutine polympo_getMPTgtElmID(mpMesh, numMPs, array) & integer(c_int), value :: numMPs type(c_ptr), value :: array end subroutine - ! + !--------------------------------------------------------------------------- !> @brief set the mp lat lon is rotational or normal !> @param mpmesh(in/out) MPMesh object @@ -159,8 +176,7 @@ subroutine polympo_setMPLatLonRotatedFlag(mpMesh, isRotateFlag) & type(c_ptr), value :: mpMesh integer(c_int), value :: isRotateFlag end subroutine - - + !--------------------------------------------------------------------------- !> @brief set the MP positions array from a host array !> @param mpmesh(in/out) MPMesh object @@ -175,6 +191,7 @@ subroutine polympo_setMPPositions(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief get the MP positions array from a polympo array !> @param mpmesh(in/out) MPMesh object @@ -190,6 +207,7 @@ subroutine polympo_getMPPositions(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief set the MP positions array from a host array !> @param mpmesh(in/out) MPMesh object @@ -204,6 +222,7 @@ subroutine polympo_setMPTgtPositions(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief get the MP positions array from a polympo array !> @param mpmesh(in/out) MPMesh object @@ -220,7 +239,6 @@ subroutine polympo_getMPTgtPositions(mpMesh, nComps, numMPs, array) & type(c_ptr), value :: array end subroutine - !--------------------------------------------------------------------------- !> @brief set the MP latitude and longtitude array from a host array !> @param mpmesh(in/out) MPMesh object @@ -236,6 +254,7 @@ subroutine polympo_setMPRotLatLon(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), intent(in), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief get the MP latitude and longtitude array from a polympo array !> @param mpmesh(in/out) MPMesh object @@ -251,6 +270,7 @@ subroutine polympo_getMPRotLatLon(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief set the MP latitude and longtitude array from a host array !> @param mpmesh(in/out) MPMesh object @@ -266,6 +286,7 @@ subroutine polympo_setMPTgtRotLatLon(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), intent(in), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief get the MP latitude and longtitude array from a polympo array !> @param mpmesh(in/out) MPMesh object @@ -282,9 +303,6 @@ subroutine polympo_getMPTgtRotLatLon(mpMesh, nComps, numMPs, array) & type(c_ptr), value :: array end subroutine - - - !--------------------------------------------------------------------------- !> @brief set the Mass MP array from a host array !> @param mpmesh(in/out) MPMesh object @@ -299,6 +317,7 @@ subroutine polympo_setMPMass(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), intent(in), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief get the Mass MP array from a polympo array !> @param mpmesh(in/out) MPMesh object @@ -313,6 +332,7 @@ subroutine polympo_getMPMass(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief set the velocity MP array from a host array !> @param mpmesh(in/out) MPMesh object @@ -327,6 +347,7 @@ subroutine polympo_setMPVel(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), intent(in), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief get the velocity MP array from a polympo array !> @param mpmesh(in/out) MPMesh object @@ -341,6 +362,17 @@ subroutine polympo_getMPVel(mpMesh, nComps, numMPs, array) & integer(c_int), value :: nComps, numMPs type(c_ptr), value :: array end subroutine + + subroutine polympo_setMPStrainRate(mpMesh) & + bind(C, NAME='polympo_setMPStrainRate_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + + !getMPStrainRate + !setMPStress + !getMPStress + !--------------------------------------------------------------------------- !> @brief Enable the setting of mesh topology (number of entities and entity adjacencies). !> By default, modifying the mesh topology without calling this function first will result @@ -352,6 +384,7 @@ subroutine polympo_startMeshFill(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief Disable the modification of mesh topology (number of entities and entity adjacencies). !> Once the topology is set, call this function to prevent unexpected/accidental modifications. @@ -362,6 +395,7 @@ subroutine polympo_endMeshFill(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief check the Mesh is valid/runable in polympo !> @param mpMesh(in/out) the MPMesh is valid/created @@ -374,6 +408,7 @@ subroutine polympo_checkMeshMaxSettings(mpMesh,maxEdges,vertexDegree) & type(c_ptr), value :: mpMesh integer(c_int), value :: maxEdges, vertexDegree end subroutine + !--------------------------------------------------------------------------- !> @brief set the Mesh type to general polygonal !> modifies mesh topology polympo_startMeshFill required @@ -384,6 +419,7 @@ subroutine polympo_setMeshTypeGeneralPoly(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief set the Mesh type to CVT polygonal !> modifies mesh topology polympo_startMeshFill required @@ -394,6 +430,7 @@ subroutine polympo_setMeshTypeCVTPoly(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief set the Mesh geometry type to planar !> modifies mesh topology polympo_startMeshFill required @@ -404,6 +441,7 @@ subroutine polympo_setMeshGeomTypePlanar(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief set the Mesh geometry type to spherical !> modifies mesh topology polympo_startMeshFill required @@ -414,6 +452,7 @@ subroutine polympo_setMeshGeomTypeSpherical(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief set the Mesh sphere radius !> modifies mesh topology polympo_startMeshFill required @@ -425,6 +464,7 @@ subroutine polympo_setMeshSphereRadius(mpMesh,sphereRadius) & type(c_ptr), value :: mpMesh real(c_double), value :: sphereRadius end subroutine + !--------------------------------------------------------------------------- !> @brief set the number of vetices of the mesh !> modifies mesh topology polympo_startMeshFill required @@ -437,6 +477,7 @@ subroutine polympo_setMeshNumVtxs(mpMesh,numVtxs) & type(c_ptr), value :: mpMesh integer(c_int), value :: numVtxs end subroutine + !--------------------------------------------------------------------------- !> @brief get the number of mesh vertices !> @param mpMesh(in) mpMesh object @@ -448,6 +489,7 @@ subroutine polympo_getMeshNumVtxs(mpMesh, numVtx) & type(c_ptr), value :: mpMesh integer(c_int), intent(inout) :: numVtx end subroutine + !--------------------------------------------------------------------------- !> @brief set the number of elements of the mesh !> modifies mesh topology polympo_startMeshFill required @@ -460,6 +502,7 @@ subroutine polympo_setMeshNumElms(mpMesh,numElms) & type(c_ptr), value :: mpMesh integer(c_int), value :: numElms end subroutine + !--------------------------------------------------------------------------- !> @brief get the number of mesh elements !> @param mpMesh(in) mpMesh object @@ -471,6 +514,7 @@ subroutine polympo_getMeshNumElms(mpMesh, numElm) & type(c_ptr), value :: mpMesh integer(c_int), intent(inout) :: numElm end subroutine + !--------------------------------------------------------------------------- !> @brief set the polympo mesh number of edges per element !> modifies mesh topology polympo_startMeshFill required @@ -485,6 +529,7 @@ subroutine polympo_setMeshNumEdgesPerElm(mpMesh, nCells, nEdgesOnCell) & integer(c_int), value :: nCells type(c_ptr), intent(in), value :: nEdgesOnCell end subroutine + !--------------------------------------------------------------------------- !> @brief set the polympo mesh element to vertices connectivity !> modifies mesh topology polympo_startMeshFill required @@ -499,6 +544,7 @@ subroutine polympo_setMeshElm2VtxConn(mpMesh, maxEdges, nCells, verticesOnCell) integer(c_int), value :: maxEdges, nCells type(c_ptr), intent(in), value :: verticesOnCell end subroutine + !--------------------------------------------------------------------------- !> @brief set the polympo mesh element to elements connectivity !> modifies mesh topology polympo_startMeshFill required @@ -513,6 +559,51 @@ subroutine polympo_setMeshElm2ElmConn(mpMesh, maxEdges, nCells, cellsOnCell) & integer(c_int), value :: maxEdges, nCells type(c_ptr), intent(in), value :: cellsOnCell end subroutine + + !--------------------------------------------------------------------------- + !> @brief set the owning process array + !> @param mpmesh(in/out) MPMesh object + !> @param nCells(in) number of cells + !> @param array(in) input mesh cell to process array + !--------------------------------------------------------------------------- + subroutine polympo_setOwningProc(mpMesh, nCells, array) & + bind(C, NAME='polympo_setOwningProc_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nCells + type(c_ptr), intent(in), value :: array + end subroutine + + subroutine polympo_setOwningProcVertex(mpMesh, nVertices, array) & + bind(C, NAME='polympo_setOwningProcVertex_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nVertices + type(c_ptr), intent(in), value :: array + end subroutine + + !--------------------------------------------------------------------------- + !> @brief set the owning process array + !> @param mpmesh(in/out) MPMesh object + !> @param nCells(in) number of cells + !> @param array(in) input mesh cell to process array + !--------------------------------------------------------------------------- + subroutine polympo_setElmGlobal(mpMesh, nCells, array) & + bind(C, NAME='polympo_setElmGlobal_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nCells + type(c_ptr), intent(in), value :: array + end subroutine + + subroutine polympo_setVtxGlobal(mpMesh, nVertices, array) & + bind(C, NAME='polympo_setVtxGlobal_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nVertices + type(c_ptr), intent(in), value :: array + end subroutine + !--------------------------------------------------------------------------- !> @brief get enum for vertex mesh fields !--------------------------------------------------------------------------- @@ -520,6 +611,7 @@ integer function polympo_getMeshFVtxType() & bind(C, NAME='polympo_getMeshFVtxType_f') use :: iso_c_binding end function + !--------------------------------------------------------------------------- !> @brief get enum for element mesh fields !--------------------------------------------------------------------------- @@ -541,6 +633,7 @@ subroutine polympo_setMeshVtxCoords(mpMesh, nVertices, xArray, yArray, zArray) & integer(c_int), value :: nVertices type(c_ptr), intent(in), value :: xArray, yArray, zArray end subroutine + !--------------------------------------------------------------------------- !> @brief get the polympo mesh vertices coordinates !> @param mpmesh(in/out) MPMesh object @@ -554,6 +647,7 @@ subroutine polympo_getMeshVtxCoords(mpMesh, nVertices, xArray, yArray, zArray) & integer(c_int), value :: nVertices type(c_ptr), value :: xArray, yArray, zArray end subroutine + !--------------------------------------------------------------------------- !> @brief set the polympo mesh vertices latitude and longitude !> @param mpmesh(in/out) MPMesh object @@ -567,6 +661,7 @@ subroutine polympo_setMeshVtxRotLat(mpMesh, nVertices, latitude) & integer(c_int), value :: nVertices type(c_ptr), intent(in), value :: latitude end subroutine + !--------------------------------------------------------------------------- !> @brief get the polympo mesh vertices latitude and longitude !> @param mpmesh(in/out) MPMesh object @@ -580,59 +675,6 @@ subroutine polympo_getMeshVtxRotLat(mpMesh, nVertices, latitude) & integer(c_int), value :: nVertices type(c_ptr), value :: latitude end subroutine - !--------------------------------------------------------------------------- - !> @brief set the polympo mesh elements/cells center - !> @param mpmesh(in/out) MPMesh object - !> @param nElements(in) length of array in, use for assertion - !> @param x/y/zArray(in) the 1D arrays of element centers coords - !--------------------------------------------------------------------------- - subroutine polympo_setMeshElmCenter(mpMesh, nCells, xArray, yArray, zArray) & - bind(C, NAME='polympo_setMeshElmCenter_f') - use :: iso_c_binding - type(c_ptr), value :: mpMesh - integer(c_int), value :: nCells - type(c_ptr), intent(in), value :: xArray, yArray, zArray - end subroutine - !--------------------------------------------------------------------------- - !> @brief get the polympo mesh elements/cells center - !> @param mpmesh(in/out) MPMesh object - !> @param nElements(in) length of array in, use for assertion - !> @param x/y/zArray(in/out) the 1D arrays of element centers coords - !--------------------------------------------------------------------------- - subroutine polympo_getMeshElmCenter(mpMesh, nCells, xArray, yArray, zArray) & - bind(C, NAME='polympo_getMeshElmCenter_f') - use :: iso_c_binding - type(c_ptr), value :: mpMesh - integer(c_int), value :: nCells - type(c_ptr), value :: xArray, yArray, zArray - end subroutine - - !--------------------------------------------------------------------------- - !> @brief set the polympo dual mesh triangle area - !> @param mpmesh(in/out) MPMesh object - !> @param nVertices(in) length of array in, use for assertion - !> @param Array(in) 1D array of area of dual triangle elements - !--------------------------------------------------------------------------- - subroutine polympo_setMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) & - bind(C, NAME='polympo_setMeshDualTriangleArea_f') - use :: iso_c_binding - type(c_ptr), value :: mpMesh - integer(c_int), value :: nVertices - type(c_ptr), intent(in), value :: areaTriangle - end subroutine - !--------------------------------------------------------------------------- - !> @brief get the polympo mesh dual mesh triangle area - !> @param mpmesh(in/out) MPMesh object - !> @param nVetices(in) length of array in, use for assertion - !> @param Array(in/out) 1D array of area of dual triangle elements - !--------------------------------------------------------------------------- - subroutine polympo_getMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) & - bind(C, NAME='polympo_getMeshDualTriangleArea_f') - use :: iso_c_binding - type(c_ptr), value :: mpMesh - integer(c_int), value :: nVertices - type(c_ptr), value :: areaTriangle - end subroutine !--------------------------------------------------------------------------- !> @brief set the vertices velocity from a host array @@ -648,6 +690,7 @@ subroutine polympo_setMeshVtxVel(mpMesh, nVertices, uVel, vVel) & integer(c_int), value :: nVertices type(c_ptr), intent(in), value :: uVel, vVel end subroutine + !--------------------------------------------------------------------------- !> @brief get the vertices velocity from polyMPO !> @param mpmesh(in/out) MPMesh object @@ -664,6 +707,7 @@ subroutine polympo_getMeshVtxVel(mpMesh, nVertices, uVel, vVel) & integer(c_int), value :: nVertices type(c_ptr), value :: uVel, vVel end subroutine + !--------------------------------------------------------------------------- !> @brief set the vertices mass from a host array !> @param mpmesh(in/out) MPMesh object @@ -677,6 +721,7 @@ subroutine polympo_setMeshVtxMass(mpMesh, nVertices, vtxMass) & integer(c_int), value :: nVertices type(c_ptr), intent(in), value :: vtxMass end subroutine + !--------------------------------------------------------------------------- !> @brief get the vertices mass from polyMPO !> @param mpmesh(in/out) MPMesh object @@ -690,6 +735,7 @@ subroutine polympo_getMeshVtxMass(mpMesh, nVertices, vtxMass) & integer(c_int), value :: nVertices type(c_ptr), value :: vtxMass end subroutine + !--------------------------------------------------------------------------- !> @brief set the mesh elements mass from a host array !> @param mpmesh(in/out) MPMesh object @@ -703,6 +749,7 @@ subroutine polympo_setMeshElmMass(mpMesh, nCells, elmMass) & integer(c_int), value :: nCells type(c_ptr), intent(in), value :: elmMass end subroutine + !--------------------------------------------------------------------------- !> @brief get the mesh element mass from polyMPO !> @param mpmesh(in/out) MPMesh object @@ -716,6 +763,7 @@ subroutine polympo_getMeshElmMass(mpMesh, nCells, elmMass) & integer(c_int), value :: nCells type(c_ptr), value :: elmMass end subroutine + !--------------------------------------------------------------------------- !> @brief set the spherical velocity increment mesh array !> from a host array @@ -731,6 +779,7 @@ subroutine polympo_setMeshVtxOnSurfVeloIncr(mpMesh, nComps, nVertices, array) & integer(c_int), value :: nComps, nVertices type(c_ptr), intent(in), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief get the spherical velocity increment mesh array !> from a polympo array @@ -747,6 +796,7 @@ subroutine polympo_getMeshVtxOnSurfVeloIncr(mpMesh, nComps, nVertices, array) & integer(c_int), value :: nComps, nVertices type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief set the spherical displacement increment mesh array !> from a host array @@ -762,6 +812,7 @@ subroutine polympo_setMeshVtxOnSurfDispIncr(mpMesh, nComps, nVertices, array) & integer(c_int), value :: nComps, nVertices type(c_ptr), intent(in), value :: array end subroutine + !--------------------------------------------------------------------------- !> @brief get the spherical displacement increment mesh array !> from a polympo array @@ -778,42 +829,77 @@ subroutine polympo_getMeshVtxOnSurfDispIncr(mpMesh, nComps, nVertices, array) & integer(c_int), value :: nComps, nVertices type(c_ptr), value :: array end subroutine + !--------------------------------------------------------------------------- - !> @brief set the owning process array + !> @brief set the polympo mesh elements/cells center !> @param mpmesh(in/out) MPMesh object - !> @param nCells(in) number of cells - !> @param array(in) input mesh cell to process array + !> @param nElements(in) length of array in, use for assertion + !> @param x/y/zArray(in) the 1D arrays of element centers coords !--------------------------------------------------------------------------- - subroutine polympo_setOwningProc(mpMesh, nCells, array) & - bind(C, NAME='polympo_setOwningProc_f') + subroutine polympo_setMeshElmCenter(mpMesh, nCells, xArray, yArray, zArray) & + bind(C, NAME='polympo_setMeshElmCenter_f') use :: iso_c_binding type(c_ptr), value :: mpMesh integer(c_int), value :: nCells - type(c_ptr), intent(in), value :: array + type(c_ptr), intent(in), value :: xArray, yArray, zArray end subroutine + !--------------------------------------------------------------------------- - !> @brief set the owning process array + !> @brief get the polympo mesh elements/cells center !> @param mpmesh(in/out) MPMesh object - !> @param nCells(in) number of cells - !> @param array(in) input mesh cell to process array + !> @param nElements(in) length of array in, use for assertion + !> @param x/y/zArray(in/out) the 1D arrays of element centers coords !--------------------------------------------------------------------------- - subroutine polympo_setElmGlobal(mpMesh, nCells, array) & - bind(C, NAME='polympo_setElmGlobal_f') + subroutine polympo_getMeshElmCenter(mpMesh, nCells, xArray, yArray, zArray) & + bind(C, NAME='polympo_getMeshElmCenter_f') use :: iso_c_binding type(c_ptr), value :: mpMesh integer(c_int), value :: nCells - type(c_ptr), intent(in), value :: array + type(c_ptr), value :: xArray, yArray, zArray end subroutine - subroutine polympo_setVtxGlobal(mpMesh, nVertices, array) & - bind(C, NAME='polympo_setVtxGlobal_f') + + !--------------------------------------------------------------------------- + !> @brief set the polympo dual mesh triangle area + !> @param mpmesh(in/out) MPMesh object + !> @param nVertices(in) length of array in, use for assertion + !> @param Array(in) 1D array of area of dual triangle elements + !--------------------------------------------------------------------------- + subroutine polympo_setMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) & + bind(C, NAME='polympo_setMeshDualTriangleArea_f') use :: iso_c_binding type(c_ptr), value :: mpMesh integer(c_int), value :: nVertices - type(c_ptr), intent(in), value :: array + type(c_ptr), intent(in), value :: areaTriangle + end subroutine + + !--------------------------------------------------------------------------- + !> @brief get the polympo mesh dual mesh triangle area + !> @param mpmesh(in/out) MPMesh object + !> @param nVetices(in) length of array in, use for assertion + !> @param Array(in/out) 1D array of area of dual triangle elements + !--------------------------------------------------------------------------- + subroutine polympo_getMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) & + bind(C, NAME='polympo_getMeshDualTriangleArea_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nVertices + type(c_ptr), value :: areaTriangle + end subroutine + + subroutine polympo_setGnomonicProjection(mpMesh) & + bind(C, NAME='polympo_setGnomonicProjection_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + + subroutine polyMPO_setTanLatVertexRotatedOverRadius(mpMesh, nVertices, array) & + bind(C, NAME=' polyMPO_setTanLatVertexRotatedOverRadius_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + integer(c_int), value :: nVertices + type(c_ptr), value :: array end subroutine - - !--------------------------------------------------------------------------- !> @brief calculate the MPs from given mesh vertices rotational latitude !> longitude, update the MP slices @@ -825,37 +911,34 @@ subroutine polympo_push(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine - - + + subroutine polympo_push_ahead(mpMesh) & + bind(C, NAME='polympo_push_ahead_f') + use :: iso_c_binding + type(c_ptr), value :: mpMesh + end subroutine + !--------------------------------------------------------------------------- !> @brief calculate the MPs from given mesh vertices rotational latitude + !--------------------------------------------------------------------------- logical function polympo_push1P(mpMesh) & bind(C, NAME='polympo_push1P_f') use :: iso_c_binding type(c_ptr), value :: mpMesh end function - subroutine polympo_push_ahead(mpMesh) & - bind(C, NAME='polympo_push_ahead_f') - use :: iso_c_binding - type(c_ptr), value :: mpMesh - end subroutine - subroutine polympo_push_swap(mpMesh) & bind(C, NAME='polympo_push_swap_f') use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine - + subroutine polympo_push_swap_pos(mpMesh) & bind(C, NAME='polympo_push_swap_pos_f') use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine - - - !--------------------------------------------------------------------------- !> @brief start the reconstruction of MP Mass to Mesh Vertices !> @param mpmesh(in/out) MPMesh object @@ -869,6 +952,7 @@ subroutine polympo_setReconstructionOfMass(mpMesh, order, meshEntType) & type(c_ptr), value :: mpMesh integer(c_int), value :: order, meshEntType end subroutine + !--------------------------------------------------------------------------- !> @brief start the reconstruction of MP Velocity to Mesh Vertices !> @param mpmesh(in/out) MPMesh object @@ -882,6 +966,7 @@ subroutine polympo_setReconstructionOfVel(mpMesh, order, meshEntType) & type(c_ptr), value :: mpMesh integer(c_int), value :: order, meshEntType end subroutine + !--------------------------------------------------------------------------- !> @brief start the reconstruction of MP Strain Rate to Mesh Vertices !> @param mpmesh(in/out) MPMesh object @@ -895,6 +980,7 @@ subroutine polympo_setReconstructionOfStrainRate(mpMesh, order, meshEntType) & type(c_ptr), value :: mpMesh integer(c_int), value :: order, meshEntType end subroutine + !--------------------------------------------------------------------------- !> @brief start the reconstruction of MP Stress to Mesh Vertices !> @param mpmesh(in/out) MPMesh object @@ -909,39 +995,24 @@ subroutine polympo_setReconstructionOfStress(mpMesh, order, meshEntType) & integer(c_int), value :: order, meshEntType end subroutine - - subroutine polympo_vtxSubAssemblyIceArea(mpMesh, vtxPerElm, nCells, comp, array) & - bind(C, NAME='polympo_vtxSubAssemblyIceArea_f') + subroutine polympo_reconstruct_coeff_with_MPI(mpMesh) & + bind(C, NAME='polympo_reconstruct_coeff_with_MPI_f') use :: iso_c_binding type(c_ptr), value :: mpMesh - integer(c_int), value :: vtxPerElm, nCells, comp - type(c_ptr), value :: array end subroutine - subroutine polympo_vtxSubAssemblyVelocity(mpMesh, vtxPerElm, nCells, comp, array) & - bind(C, NAME='polympo_vtxSubAssemblyVelocity_f') + subroutine polympo_reconstruct_iceArea_with_MPI(mpMesh) & + bind(C, NAME='polympo_reconstruct_iceArea_with_MPI_f') use :: iso_c_binding type(c_ptr), value :: mpMesh - integer(c_int), value :: vtxPerElm, nCells, comp - type(c_ptr), value :: array end subroutine - - subroutine polympo_subAssemblyCoeffs(mpMesh, vtxPerElm, nCells, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44) & - bind(C, NAME='polympo_subAssemblyCoeffs_f') - use :: iso_c_binding - type(c_ptr), value :: mpMesh - integer(c_int), value :: vtxPerElm, nCells - type(c_ptr), value :: m11, m12, m13, m14, m22, m23, m24, m33, m34, m44 - end subroutine - - subroutine polympo_regularize_and_solve_matrix(mpMesh, nVertices, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44) & - bind(C, NAME='polympo_regularize_and_solve_matrix_f') + + subroutine polympo_reconstruct_velocity_with_MPI(mpMesh) & + bind(C, NAME='polympo_reconstruct_velocity_with_MPI_f') use :: iso_c_binding type(c_ptr), value :: mpMesh - integer(c_int), value :: nVertices - type(c_ptr), value :: m11, m12, m13, m14, m22, m23, m24, m33, m34, m44 end subroutine - + !--------------------------------------------------------------------------- !> @brief directly call the reconstruct of the MP fields to mesh fields !> @param mpmesh(in/out) MPMesh object @@ -951,18 +1022,21 @@ subroutine polympo_applyReconstruction(mpMesh) & use :: iso_c_binding type(c_ptr), value :: mpMesh end subroutine + !--------------------------------------------------------------------------- !> @brief enable timing functions !--------------------------------------------------------------------------- subroutine polympo_enableTiming() bind(C, NAME='polympo_enableTiming_f') use :: iso_c_binding end subroutine + !--------------------------------------------------------------------------- !> @brief print summary of timings collected !--------------------------------------------------------------------------- subroutine polympo_summarizeTime() bind(C, NAME='polympo_summarizeTime_f') use :: iso_c_binding end subroutine + !--------------------------------------------------------------------------- !> @brief set verbosity of timings collected !--------------------------------------------------------------------------- @@ -971,12 +1045,6 @@ subroutine polympo_setTimingVerbosity(v) bind(C, NAME='polympo_setTimingVerbosit integer(c_int), value :: v end subroutine - integer function polympo_getMPCount(mpMesh) bind(C, name="polympo_getMPCount_f") - use :: iso_c_binding - implicit none - type(c_ptr), value :: mpMesh - end function - end interface contains !--------------------------------------------------------------------------- diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp index f8391d1..ef46297 100644 --- a/src/pmpo_materialPoints.hpp +++ b/src/pmpo_materialPoints.hpp @@ -125,7 +125,6 @@ class MaterialPoints { PS* MPs; int elmIDoffset = -1; int maxAppID = -1; - bool isRotatedFlag = false; Operating_Mode operating_mode; RebuildHelper rebuildFields; IntFunc getAppID; @@ -211,46 +210,47 @@ class MaterialPoints { ps::parallel_for(MPs, swap, "swap"); } void updateMPSliceAll(){ - updateMPElmID(); - updateMPSlice(); - updateMPSlice(); + updateMPElmID(); + updateMPSlice(); + updateMPSlice(); + } + + void updateRotLatLonAndXYZ2Tgt(const double radius, const bool isRotated){ + Kokkos::Timer timer; + auto curPosRotLatLon = MPs->get(); + auto tgtPosRotLatLon = MPs->get(); + auto tgtPosXYZ = MPs->get(); + auto rotLatLonIncr = MPs->get(); + //Velocity + auto velMPs = MPs->get(); + auto velIncr = MPs->get(); + + auto mpAppID = MPs->get(); + + 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 + tgtPosRotLatLon(mp,0) = rotLat; + tgtPosRotLatLon(mp,1) = rotLon; + auto geoLat = rotLat; + auto geoLon = rotLon; + if(isRotated){ + 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); + velMPs(mp,0) = velMPs(mp,0) + velIncr(mp,0); + velMPs(mp,1) = velMPs(mp,1) + velIncr(mp,1); + } + }; + ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude"); + pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds()); } - void updateRotLatLonAndXYZ2Tgt(const double radius){ - Kokkos::Timer timer; - auto curPosRotLatLon = MPs->get(); - auto tgtPosRotLatLon = MPs->get(); - auto tgtPosXYZ = MPs->get(); - auto rotLatLonIncr = MPs->get(); - //Velocity - auto velMPs = MPs->get(); - auto velIncr = MPs->get(); - - 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 - tgtPosRotLatLon(mp,0) = rotLat; - tgtPosRotLatLon(mp,1) = rotLon; - auto geoLat = rotLat; - auto geoLon = rotLon; - 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); - - velMPs(mp,0) = velMPs(mp,0) + velIncr(mp,0); - velMPs(mp,1) = velMPs(mp,1) + velIncr(mp,1); - } - }; - ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude"); - pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds()); - } template auto getData() { @@ -281,13 +281,6 @@ class MaterialPoints { PMT_ALWAYS_ASSERT(maxAppID != -1); return maxAppID; } - bool getRotatedFlag() { - return isRotatedFlag; - } - void setRotatedFlag(bool flagSet) { - isRotatedFlag = flagSet; - } - // MUTATOR template void fillData(double value);//use PS_LAMBDA fill up to 1 };// End MaterialPoints diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index 57b838f..0e3760a 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -1,87 +1,128 @@ #include "pmpo_mesh.hpp" namespace polyMPO{ - bool Mesh::checkMeshType(int meshType){ - return (meshType >mesh_unrecognized_lower && - meshType mesh_unrecognized_lower && - geomType (vtxCoordsMapEntry.second,numVtxs_); - - auto vtxRotLatMapEntry = meshFields2TypeAndString.at(MeshF_VtxRotLat); - PMT_ALWAYS_ASSERT(vtxRotLatMapEntry.first == MeshFType_VtxBased); - vtxRotLat_ = MeshFView(vtxRotLatMapEntry.second,numVtxs_); - - auto vtxVelMapEntry = meshFields2TypeAndString.at(MeshF_Vel); - PMT_ALWAYS_ASSERT(vtxVelMapEntry.first == MeshFType_VtxBased); - vtxVel_ = MeshFView(vtxVelMapEntry.second,numVtxs_); - - auto vtxMassMapEntry = meshFields2TypeAndString.at(MeshF_VtxMass); - PMT_ALWAYS_ASSERT(vtxMassMapEntry.first == MeshFType_VtxBased); - vtxMass_ = MeshFView(vtxMassMapEntry.second,numVtxs_); - - auto vtxOnSurfVeloIncrMapEntry = meshFields2TypeAndString.at(MeshF_OnSurfVeloIncr); - PMT_ALWAYS_ASSERT(vtxOnSurfVeloIncrMapEntry.first == MeshFType_VtxBased); - vtxOnSurfVeloIncr_ = MeshFView(vtxOnSurfVeloIncrMapEntry.second,numVtxs_); - - auto vtxOnSurfDispIncrMapEntry = meshFields2TypeAndString.at(MeshF_OnSurfDispIncr); - PMT_ALWAYS_ASSERT(vtxOnSurfDispIncrMapEntry.first == MeshFType_VtxBased); - vtxOnSurfDispIncr_ = MeshFView(vtxOnSurfDispIncrMapEntry.second,numVtxs_); - - auto vtxRotLatLonIncrMapEntry = meshFields2TypeAndString.at(MeshF_RotLatLonIncr); - PMT_ALWAYS_ASSERT(vtxRotLatLonIncrMapEntry.first == MeshFType_VtxBased); - vtxRotLatLonIncr_ = MeshFView(vtxRotLatLonIncrMapEntry.second,numVtxs_); - - auto dualTriangleAreaEntry = meshFields2TypeAndString.at(MeshF_DualTriangleArea); - PMT_ALWAYS_ASSERT(dualTriangleAreaEntry.first == MeshFType_VtxBased); - dualTriangleArea_ = MeshFView(dualTriangleAreaEntry.second,numVtxs_); - } - - void Mesh::setMeshElmBasedFieldSize(){ - PMT_ALWAYS_ASSERT(meshEdit_); - - auto elmMassMapEntry = meshFields2TypeAndString.at(MeshF_ElmMass); - PMT_ALWAYS_ASSERT(elmMassMapEntry.first == MeshFType_ElmBased); - elmMass_ = MeshFView(elmMassMapEntry.second,numElms_); - elmCenterXYZ_ = MeshFView(meshFields2TypeAndString.at(MeshF_ElmCenterXYZ).second,numElms_); - } - - void Mesh::computeRotLatLonIncr(){ - Kokkos::Timer timer; - PMT_ALWAYS_ASSERT(geomType_ == geom_spherical_surf); - - auto dispIncr = getMeshField(); - auto rotLatLonIncr = getMeshField(); - auto lat = getMeshField(); - auto sphereRadius = getSphereRadius(); - PMT_ALWAYS_ASSERT(sphereRadius > 0); - Kokkos::parallel_for("set nEdgesPerElm", numVtxs_, KOKKOS_LAMBDA(const int iVtx){ - // Lat [iVtx,0] = dispIncrY [iVtx,1] /R - // Lon [iVtx,1] = dispIncrX [iVtx,0] /(R*cos(lat)) - rotLatLonIncr(iVtx, 0) = dispIncr(iVtx, 1)/sphereRadius; - rotLatLonIncr(iVtx, 1) = dispIncr(iVtx, 0)/(sphereRadius * std::cos(lat(iVtx))); - }); - pumipic::RecordTime("PolyMPO_computeRotLatLonIncr", timer.seconds()); - } - - IntView Mesh::getElm2Process() { - return owningProc_; - } - - IntView Mesh::getElmGlobal() { - return globalElm_; - } - - + bool Mesh::checkMeshType(int meshType){ + return (meshType >mesh_unrecognized_lower && meshType mesh_unrecognized_lower && geomType (vtxCoordsMapEntry.second,numVtxs_); + + auto vtxRotLatMapEntry = meshFields2TypeAndString.at(MeshF_VtxRotLat); + PMT_ALWAYS_ASSERT(vtxRotLatMapEntry.first == MeshFType_VtxBased); + vtxRotLat_ = MeshFView(vtxRotLatMapEntry.second,numVtxs_); + + auto vtxVelMapEntry = meshFields2TypeAndString.at(MeshF_Vel); + PMT_ALWAYS_ASSERT(vtxVelMapEntry.first == MeshFType_VtxBased); + vtxVel_ = MeshFView(vtxVelMapEntry.second,numVtxs_); + + auto vtxMassMapEntry = meshFields2TypeAndString.at(MeshF_VtxMass); + PMT_ALWAYS_ASSERT(vtxMassMapEntry.first == MeshFType_VtxBased); + vtxMass_ = MeshFView(vtxMassMapEntry.second,numVtxs_); + + auto vtxOnSurfVeloIncrMapEntry = meshFields2TypeAndString.at(MeshF_OnSurfVeloIncr); + PMT_ALWAYS_ASSERT(vtxOnSurfVeloIncrMapEntry.first == MeshFType_VtxBased); + vtxOnSurfVeloIncr_ = MeshFView(vtxOnSurfVeloIncrMapEntry.second,numVtxs_); + + auto vtxOnSurfDispIncrMapEntry = meshFields2TypeAndString.at(MeshF_OnSurfDispIncr); + PMT_ALWAYS_ASSERT(vtxOnSurfDispIncrMapEntry.first == MeshFType_VtxBased); + vtxOnSurfDispIncr_ = MeshFView(vtxOnSurfDispIncrMapEntry.second,numVtxs_); + + auto vtxRotLatLonIncrMapEntry = meshFields2TypeAndString.at(MeshF_RotLatLonIncr); + PMT_ALWAYS_ASSERT(vtxRotLatLonIncrMapEntry.first == MeshFType_VtxBased); + vtxRotLatLonIncr_ = MeshFView(vtxRotLatLonIncrMapEntry.second,numVtxs_); + + auto dualTriangleAreaEntry = meshFields2TypeAndString.at(MeshF_DualTriangleArea); + PMT_ALWAYS_ASSERT(dualTriangleAreaEntry.first == MeshFType_VtxBased); + dualTriangleArea_ = MeshFView(dualTriangleAreaEntry.second,numVtxs_); + + auto tanLatVertexRotOverRadiusEntry = meshFields2TypeAndString.at(MeshF_TanLatVertexRotatedOverRadius); + PMT_ALWAYS_ASSERT(tanLatVertexRotOverRadiusEntry.first == MeshFType_VtxBased); + tanLatVertexRotatedOverRadius_ = MeshFView(tanLatVertexRotOverRadiusEntry.second, numVtxs_); + + } + + void Mesh::setMeshElmBasedFieldSize(){ + PMT_ALWAYS_ASSERT(meshEdit_); + + auto elmMassMapEntry = meshFields2TypeAndString.at(MeshF_ElmMass); + PMT_ALWAYS_ASSERT(elmMassMapEntry.first == MeshFType_ElmBased); + elmMass_ = MeshFView(elmMassMapEntry.second,numElms_); + + elmCenterXYZ_ = MeshFView(meshFields2TypeAndString.at(MeshF_ElmCenterXYZ).second, numElms_); + + elmCenterGnomProj_= MeshFView(meshFields2TypeAndString.at(MeshF_ElmCenterGnomProj).second, numElms_); + + vtxGnomProj_ = MeshFView(meshFields2TypeAndString.at(MeshF_VtxGnomProj).second, numElms_); + } + + void Mesh::computeRotLatLonIncr(){ + Kokkos::Timer timer; + PMT_ALWAYS_ASSERT(geomType_ == geom_spherical_surf); + + auto dispIncr = getMeshField(); + auto rotLatLonIncr = getMeshField(); + auto lat = getMeshField(); + auto sphereRadius = getSphereRadius(); + PMT_ALWAYS_ASSERT(sphereRadius > 0); + Kokkos::parallel_for("set nEdgesPerElm", numVtxs_, KOKKOS_LAMBDA(const int iVtx){ + // Lat [iVtx,0] = dispIncrY [iVtx,1] /R + // Lon [iVtx,1] = dispIncrX [iVtx,0] /(R*cos(lat)) + rotLatLonIncr(iVtx, 0) = dispIncr(iVtx, 1)/sphereRadius; + rotLatLonIncr(iVtx, 1) = dispIncr(iVtx, 0)/(sphereRadius * std::cos(lat(iVtx))); + }); + pumipic::RecordTime("PolyMPO_computeRotLatLonIncr", timer.seconds()); + } + + void Mesh::setGnomonicProjection(bool isRotated){ + std::cout<<__FUNCTION__<(); + auto gnomProjElmCenter = getMeshField(); + + auto vtxCoords = getMeshField(); + auto elmCenters = getMeshField(); + auto elm2VtxConn = getElm2VtxConn(); + + Kokkos::parallel_for("setGnomprojCenter", numElms_, KOKKOS_LAMBDA(const int iElm){ + Vec3d elmCenter(elmCenters(iElm, 0), elmCenters(iElm, 1), elmCenters(iElm, 2)); + if(isRotated){ + elmCenter[0] = - elmCenters(iElm, 2); + elmCenter[2] = elmCenters(iElm, 0); + } + auto cos2LatR = elmCenter[0]*elmCenter[0] + elmCenter[1]*elmCenter[1]; + auto invR = 1.0/ sqrt(cos2LatR + elmCenter[2]*elmCenter[2]); + auto cosLatR = sqrt(cos2LatR); + + gnomProjElmCenter(iElm, 0) = elmCenter[1]/cosLatR; + gnomProjElmCenter(iElm, 1) = elmCenter[0]/cosLatR; + gnomProjElmCenter(iElm, 2) = invR*elmCenter[2]; + gnomProjElmCenter(iElm, 3) = invR*cosLatR; + + int nVtxE = elm2VtxConn(iElm,0); + for(int i=0; i struct meshFieldToType < MeshF_ElmMass > { using type = Ko template <> struct meshFieldToType < MeshF_OnSurfVeloIncr > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_OnSurfDispIncr > { using type = Kokkos::View; }; template <> struct meshFieldToType < MeshF_RotLatLonIncr > { using type = Kokkos::View; }; +template <> struct meshFieldToType < MeshF_VtxGnomProj > { using type = Kokkos::View; }; +template <> struct meshFieldToType < MeshF_ElmCenterGnomProj > { using type = Kokkos::View; }; +template <> struct meshFieldToType < MeshF_TanLatVertexRotatedOverRadius > { using type = Kokkos::View; }; template using MeshFView = typename meshFieldToType::type; -const std::map> meshFields2TypeAndString = - {{MeshF_Invalid, {MeshFType_Invalid,"MeshField_InValid!"}}, - {MeshF_Unsupported, {MeshFType_Unsupported,"MeshField_Unsupported"}}, - {MeshF_VtxCoords, {MeshFType_VtxBased,"MeshField_VerticesCoords"}}, - {MeshF_VtxRotLat, {MeshFType_VtxBased,"MeshField_VerticesLatitude"}}, - {MeshF_ElmCenterXYZ, {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}}, - {MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}}, - {MeshF_Vel, {MeshFType_VtxBased,"MeshField_Velocity"}}, - {MeshF_VtxMass, {MeshFType_VtxBased,"MeshField_VerticesMass"}}, - {MeshF_ElmMass, {MeshFType_ElmBased,"MeshField_ElementsMass"}}, - {MeshF_OnSurfVeloIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceVelocityIncrement"}}, - {MeshF_OnSurfDispIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}}, - {MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}}}; +const std::map> meshFields2TypeAndString = { + {MeshF_Invalid, {MeshFType_Invalid,"MeshField_InValid!"}}, + {MeshF_Unsupported, {MeshFType_Unsupported,"MeshField_Unsupported"}}, + {MeshF_VtxCoords, {MeshFType_VtxBased,"MeshField_VerticesCoords"}}, + {MeshF_VtxRotLat, {MeshFType_VtxBased,"MeshField_VerticesLatitude"}}, + {MeshF_ElmCenterXYZ, {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}}, + {MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}}, + {MeshF_Vel, {MeshFType_VtxBased,"MeshField_Velocity"}}, + {MeshF_VtxMass, {MeshFType_VtxBased,"MeshField_VerticesMass"}}, + {MeshF_ElmMass, {MeshFType_ElmBased,"MeshField_ElementsMass"}}, + {MeshF_OnSurfVeloIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceVelocityIncrement"}}, + {MeshF_OnSurfDispIncr, {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}}, + {MeshF_RotLatLonIncr, {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}}, + {MeshF_VtxGnomProj, {MeshFType_ElmBased,"MeshField_VertexGnomonicProjection"}}, + {MeshF_ElmCenterGnomProj,{MeshFType_ElmBased,"MeshField_ElementCenterGnomonicprojection"}}, + {MeshF_TanLatVertexRotatedOverRadius, {MeshFType_VtxBased,"MeshField_TanLatVertexRotatedOverRadius"}}, +}; enum mesh_type {mesh_unrecognized_lower = -1, mesh_general_polygonal, //other meshes @@ -88,6 +97,7 @@ class Mesh { IntVtx2ElmView elm2VtxConn_; IntElm2ElmView elm2ElmConn_; IntView owningProc_; + IntView owningProcVertex_; IntView globalElm_; IntView globalVtx_; //start of meshFields @@ -101,7 +111,12 @@ class Mesh { MeshFView vtxOnSurfVeloIncr_; MeshFView vtxOnSurfDispIncr_; MeshFView vtxRotLatLonIncr_; + //GnomonicProjection + MeshFView vtxGnomProj_; + MeshFView elmCenterGnomProj_; //DoubleMat2DView vtxStress_; + MeshFView tanLatVertexRotatedOverRadius_; + bool isRotatedFlag = false; public: Mesh(){}; @@ -125,13 +140,21 @@ class Mesh { setMeshElmBasedFieldSize(); meshEdit_ = false; vtxCoords_ = vtxCoords; - } + } bool meshEditable(){ return meshEdit_; } bool checkMeshType(int meshType); bool checkGeomType(int geomType); - IntView getElm2Process(); + void setOwningProc(IntView owningProc){ + PMT_ALWAYS_ASSERT(meshEdit_); + owningProc_ = owningProc; + } + void setOwningProcVertex(IntView owningProcVertex){ + owningProcVertex_ = owningProcVertex; + } + IntView getElm2Process() {return owningProc_;} + IntView getVtx2Process() {return owningProcVertex_;} mesh_type getMeshType() { return meshType_; } geom_type getGeomType() { return geomType_; } @@ -161,15 +184,23 @@ class Mesh { elm2VtxConn_ = elm2VtxConn; } void setElm2ElmConn(IntElm2ElmView elm2ElmConn) {PMT_ALWAYS_ASSERT(meshEdit_); elm2ElmConn_ = elm2ElmConn; } - void setOwningProc(IntView owningProc) {PMT_ALWAYS_ASSERT(meshEdit_); - owningProc_ = owningProc; } - + + void setElmGlobal(IntView globalElm) {globalElm_ = globalElm;} - IntView getElmGlobal(); void setVtxGlobal(IntView globalVtx) {globalVtx_ = globalVtx;} + IntView getElmGlobal() {return globalElm_;} IntView getVtxGlobal() {return globalVtx_;} + void setGnomonicProjection(bool isRotated); + void computeRotLatLonIncr(); + + bool getRotatedFlag() { + return isRotatedFlag; + } + void setRotatedFlag(bool flagSet) { + isRotatedFlag = flagSet; + } }; template @@ -212,6 +243,15 @@ auto Mesh::getMeshField(){ else if constexpr (index==MeshF_RotLatLonIncr){ return vtxRotLatLonIncr_; } + else if constexpr (index==MeshF_VtxGnomProj){ + return vtxGnomProj_; + } + else if constexpr (index==MeshF_ElmCenterGnomProj){ + return elmCenterGnomProj_; + } + else if constexpr (index==MeshF_TanLatVertexRotatedOverRadius){ + return tanLatVertexRotatedOverRadius_; + } fprintf(stderr,"Mesh Field Index error!\n"); exit(1); } @@ -225,6 +265,18 @@ void Mesh::fillMeshField(int size, int numEntries, double val){ }); } +KOKKOS_INLINE_FUNCTION +void computeGnomonicProjectionAtPoint(const Vec3d& Coord, + const Kokkos::View>& gnomProjElmCenter_sub, + double& outX, double& outY){ + const double iDen = 1.0 / (gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(3) * Coord[0] + + gnomProjElmCenter_sub(0) * gnomProjElmCenter_sub(3) * Coord[1] + + gnomProjElmCenter_sub(2) * Coord[2]); + outX = iDen * (Coord[1] * gnomProjElmCenter_sub(1) - + Coord[0] * gnomProjElmCenter_sub(0)); + outY = iDen * (Coord[2] * gnomProjElmCenter_sub(3) - Coord[1] * gnomProjElmCenter_sub(2) * gnomProjElmCenter_sub(0) - + Coord[0] * gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(2)); +} } #endif diff --git a/src/pmpo_utils.hpp b/src/pmpo_utils.hpp index 76c2a7d..9f447d7 100644 --- a/src/pmpo_utils.hpp +++ b/src/pmpo_utils.hpp @@ -51,7 +51,7 @@ class Vec2d { //constructors KOKKOS_INLINE_FUNCTION Vec2d():coords_{0.0, 0.0}{} - + ~Vec2d() = default; KOKKOS_INLINE_FUNCTION @@ -63,31 +63,31 @@ class Vec2d { //operators KOKKOS_INLINE_FUNCTION double operator[](int i) const { return coords_[i]; } - + KOKKOS_INLINE_FUNCTION double &operator[](int i) { return coords_[i]; } KOKKOS_INLINE_FUNCTION Vec2d operator-() { return Vec2d(-coords_[0], -coords_[1]); } - + KOKKOS_INLINE_FUNCTION Vec2d operator+(Vec2d v) { return Vec2d(coords_[0] + v[0], coords_[1] + v[1]); } - + KOKKOS_INLINE_FUNCTION Vec2d operator-(Vec2d v) { return Vec2d(coords_[0] - v[0], coords_[1] - v[1]); } - + KOKKOS_INLINE_FUNCTION Vec2d operator*(double scalar) const { return Vec2d(coords_[0]*scalar, coords_[1]*scalar); } - + KOKKOS_INLINE_FUNCTION double dot(Vec2d v) { return coords_[0]*v[0]+coords_[1]*v[1]; } - + KOKKOS_INLINE_FUNCTION double cross(Vec2d v) { return (coords_[0]*v[1] - coords_[1]*v[0]); } - + KOKKOS_INLINE_FUNCTION double magnitude() {return sqrt(coords_[0]*coords_[0] + coords_[1]*coords_[1]); } @@ -113,10 +113,10 @@ class Vec3d{ //operators KOKKOS_INLINE_FUNCTION double operator[](int i) const { return coords_[i]; } - + KOKKOS_INLINE_FUNCTION double &operator[](int i) { return coords_[i]; } - + KOKKOS_INLINE_FUNCTION Vec3d operator+(const Vec3d& v) const { return Vec3d(coords_[0] + v.coords_[0], coords_[1] + v.coords_[1], coords_[2] + v.coords_[2]); @@ -126,7 +126,7 @@ class Vec3d{ Vec3d operator-(const Vec3d& v) const { return Vec3d(coords_[0] - v.coords_[0], coords_[1] - v.coords_[1], coords_[2] - v.coords_[2]); } - + KOKKOS_INLINE_FUNCTION Vec3d operator-() { return Vec3d(-coords_[0], -coords_[1], -coords_[2]); } @@ -160,7 +160,7 @@ class Vec4d{ vec4d_t coords_; public: - + //constructors KOKKOS_INLINE_FUNCTION Vec4d():coords_{0.0, 0.0, 0.0, 0.0}{} @@ -184,43 +184,126 @@ class Vec4d{ return Vec4d(coords_[0] + v.coords_[0], coords_[1] + v.coords_[1], coords_[2] + v.coords_[2], coords_[3] + v.coords_[3]); } - + KOKKOS_INLINE_FUNCTION Vec4d operator-(const Vec4d& v) const { return Vec4d(coords_[0] - v.coords_[0], coords_[1] - v.coords_[1], coords_[2] - v.coords_[2], coords_[3] - v.coords_[3]); - } - + } + KOKKOS_INLINE_FUNCTION Vec4d operator-() { return Vec4d(-coords_[0], -coords_[1], -coords_[2], -coords_[3]); } - + KOKKOS_INLINE_FUNCTION Vec4d operator*(double scalar) const { return Vec4d(coords_[0] * scalar, coords_[1] * scalar, coords_[2] * scalar, coords_[3] * scalar); } - + KOKKOS_INLINE_FUNCTION double dot(const Vec4d& v) const { return coords_[0] * v.coords_[0] + coords_[1] * v.coords_[1] + coords_[2] * v.coords_[2] + coords_[3] * v.coords_[3]; } - + KOKKOS_INLINE_FUNCTION double magnitude() const { return std::sqrt(coords_[0] * coords_[0] + coords_[1] * coords_[1] + coords_[2] * coords_[2] + coords_[3] * coords_[3]); } }; -class Matrix4d { - +class Matrix3d { + private: - + double data_[3][3]; + + public: + + // Default constructor: zero initialize + KOKKOS_INLINE_FUNCTION + Matrix3d() { + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + data_[i][j] = 0.0; + } + + // Constructor from 3 Vec3d rows + KOKKOS_INLINE_FUNCTION + Matrix3d(Vec3d v0, Vec3d v1, Vec3d v2) { + for (int i = 0; i < 3; ++i) { + data_[0][i] = v0[i]; + data_[1][i] = v1[i]; + data_[2][i] = v2[i]; + } + } + // Element access + KOKKOS_INLINE_FUNCTION + double& operator()(int i, int j) { return data_[i][j]; } + + KOKKOS_INLINE_FUNCTION + const double& operator()(int i, int j) const { return data_[i][j]; } + + //Matrix multiplication + KOKKOS_INLINE_FUNCTION + Matrix3d operator*(const Matrix3d& B)const { + Matrix3d C; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + double sum = 0.0; + for (int k = 0; k < 3; ++k) { + sum += data_[i][k] * B(k, j); + } + C(i, j) = sum; + } + } + return C; + } + + // Vector Matrix multiplication: y = x_Transpose*A + KOKKOS_INLINE_FUNCTION + Vec3d operator*(const Vec3d& x) const { + Vec3d y; + for (int j = 0; j < 3; ++j) { + double sum = 0.0; + for (int i = 0; i < 3; ++i) { + sum += x[i] * data_[i][j]; + } + y[j] = sum; + } + return y; + } + + // Matrix Vector multiplication A*x + KOKKOS_INLINE_FUNCTION + Vec3d rightMultiply(const Vec3d& x) const { + Vec3d y; + for (int i = 0; i < 3; ++i) { + double sum = 0.0; + for (int j = 0; j < 3; ++j) { + sum += data_[i][j]*x[j]; + } + y[i] = sum; + } + return y; + } + + KOKKOS_INLINE_FUNCTION + Matrix3d transpose() const { + Matrix3d result; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + result(i, j) = data_[j][i]; + return result; + } +}; + +class Matrix4d { + + private: double data_[4][4]; - public: - + KOKKOS_INLINE_FUNCTION Matrix4d(Vec4d v0, Vec4d v1, Vec4d v2, Vec4d v3){ - + for (int i=0; i<4; i++){ data_[0][i] = v0[i]; data_[1][i] = v1[i]; @@ -228,7 +311,7 @@ class Matrix4d { data_[3][i] = v3[i]; } } - + //Retrieval KOKKOS_INLINE_FUNCTION double& operator()(int i, int j) { return data_[i][j];} @@ -266,7 +349,7 @@ class Matrix4d { } return traceSum; } - + //Function to regularize the matrix KOKKOS_INLINE_FUNCTION void addToDiag(double eps) { @@ -286,7 +369,7 @@ class Matrix4d { data_[i][0] *= factor; } } - + }; @@ -467,7 +550,7 @@ 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 diff --git a/src/pmpo_wachspressBasis.hpp b/src/pmpo_wachspressBasis.hpp index 1dc262b..f04c1f3 100644 --- a/src/pmpo_wachspressBasis.hpp +++ b/src/pmpo_wachspressBasis.hpp @@ -6,6 +6,126 @@ namespace polyMPO{ +// spherical interpolation of values from mesh vertices to MPs +template +void sphericalInterpolation(MPMesh& mpMesh){ + Kokkos::Timer timer; + + auto p_mesh = mpMesh.p_mesh; + auto vtxCoords = p_mesh->getMeshField(); + int numVtxs = p_mesh->getNumVertices(); + auto elm2VtxConn = p_mesh->getElm2VtxConn(); + double radius = p_mesh->getSphereRadius(); + PMT_ALWAYS_ASSERT(radius >0); + + auto p_MPs = mpMesh.p_MPs; + auto MPsPosition = p_MPs->getPositions(); + auto MPsBasis = p_MPs->getData(); + + constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice; + auto mpField = p_MPs->getData(); + + const int numEntries = mpSliceToNumEntries(); + auto meshField = p_mesh->getMeshField(); + + auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { + if(mask) { //if material point is 'active'/'enabled' + int numVtx = elm2VtxConn(elm,0); + for(int entry=0; entryparallel_for(interpolation, "interpolation"); + pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds()); +} + +KOKKOS_INLINE_FUNCTION +void wachpress_weights_grads_2D(int numVtx, const Kokkos::View>& gnom_vtx_subview, + double mpProjX, double mpProjY, double radius, double* basis, double* grad_basis){ + + double vertCoords[2][maxVtxsPerElm + 1]; + for (int i = 0; i < numVtx; ++i) { + vertCoords[0][i] = gnom_vtx_subview(i, 0); + vertCoords[1][i] = gnom_vtx_subview(i, 1); + } + vertCoords[0][numVtx] = vertCoords[0][0]; + vertCoords[1][numVtx] = vertCoords[1][0]; + + //Compute areaV and areaXV + double areaV[maxVtxsPerElm]; + double areaXV[maxVtxsPerElm]; + double xy[2] = { mpProjX, mpProjY }; + + //Helper lambda for 2D triangle area + auto triArea = [&](const double p1[2], const double p2[2], const double p3[2]) -> double { + return 0.5 * (p1[0] * (p2[1] - p3[1]) - p2[0] * (p1[1] - p3[1]) + p3[0] * (p1[1] - p2[1])); + }; + + //Special case + double p1[2] = { vertCoords[0][numVtx - 1], vertCoords[1][numVtx - 1] }; + double p2[2] = { vertCoords[0][0], vertCoords[1][0] }; + double p3[2] = { vertCoords[0][1], vertCoords[1][1] }; + areaV[0] = triArea(p1, p2, p3); + areaXV[0] = triArea(p2, xy, p3); + + for (int i = 1; i < numVtx; ++i) { + double p1[2] = { vertCoords[0][i - 1], vertCoords[1][i - 1] }; + double p2[2] = { vertCoords[0][i], vertCoords[1][i] }; + double p3[2] = { vertCoords[0][i + 1], vertCoords[1][i + 1] }; + areaV[i] = triArea(p1, p2, p3); + areaXV[i] = triArea(p2, xy, p3); + } + + double denominator = 0.0; + double derivative_sum[2] = {0.0}; + double derivative[2][maxVtxsPerElm]; + double W[maxVtxsPerElm]; + + for (int i = 0; i < numVtx; ++i){ + double product = areaV[i]; + double product_sum[2] = {0.0}; + + for (int j = 0; j < numVtx - 2; ++j) { + int ind1 = (i + j + 1) % numVtx; + product *= areaXV[ind1]; + double product_dx[2] = {areaV[i], areaV[i]}; + + for (int k = 0; k < numVtx - 2; k++){ + if (k == j) continue; + int ind2 = (i + k + 1) % numVtx; + product_dx[0] = product_dx[0] * areaXV[ind2]; + product_dx[1] = product_dx[1] * areaXV[ind2]; + } + product_dx[0] = product_dx[0] * 0.5 * (vertCoords[1][ind1+1]- vertCoords[1][ind1]); + product_dx[1] = -product_dx[1] * 0.5 * (vertCoords[0][ind1+1]- vertCoords[0][ind1]); + + product_sum[0] += product_dx[0]; + product_sum[1] += product_dx[1]; + } + W[i] = product; + denominator += product; + + derivative[0][i] = product_sum[0]; + derivative[1][i] = product_sum[1]; + + derivative_sum[0] += product_sum[0]; + derivative_sum[1] += product_sum[1]; + } + + for (int i = 0; i < numVtx; ++i){ + grad_basis[i*2 + 0] = derivative[0][i] / denominator - (W[i] / (denominator * denominator)) * derivative_sum[0]; + grad_basis[i*2 + 0] = grad_basis[i*2 + 0] / radius; + grad_basis[i*2 + 1] = derivative[1][i] / denominator - (W[i] / (denominator * denominator)) * derivative_sum[1]; + grad_basis[i*2 + 1] = grad_basis[i*2 + 1] / radius; + basis[i] = W[i] / denominator; + } +} + /** \brief calculate the basis and gradient of Basis for a give MP with its element Vtxs * * \details based on the 4.1 section from: @@ -272,7 +392,6 @@ void getBasisByAreaGblFormSpherical2(Vec3d MP, int numVtxs, Vec3d* v, calcBasis(numVtxs, a, c, basis); } - /* KOKKOS_INLINE_FUNCTION void getBasisByAreaGblForm_1(Vec2d MP, int numVtxs, Vec2d* vtxCoords, double* basis) { @@ -303,129 +422,5 @@ void getBasisByAreaGblForm_1(Vec2d MP, int numVtxs, Vec2d* vtxCoords, double* ba } */ -// spherical interpolation of values from mesh vertices to MPsi -template -void sphericalInterpolation(MPMesh& mpMesh){ - Kokkos::Timer timer; - auto p_mesh = mpMesh.p_mesh; - auto vtxCoords = p_mesh->getMeshField(); - int numVtxs = p_mesh->getNumVertices(); - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - - auto p_MPs = mpMesh.p_MPs; - auto MPsPosition = p_MPs->getPositions(); - double radius = p_mesh->getSphereRadius(); - PMT_ALWAYS_ASSERT(radius >0); - constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice; - auto mpField = p_MPs->getData(); - - const int numEntries = mpSliceToNumEntries(); - //check field correspondence - auto meshField = p_mesh->getMeshField(); - - auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { - if(mask) { //if material point is 'active'/'enabled' - Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2)); - // formating - Vec3d v3d[maxVtxsPerElm+1]; - int numVtx = elm2VtxConn(elm,0); - for(int i = 1; i<=numVtx; i++){ - v3d[i-1][0] = vtxCoords(elm2VtxConn(elm,i)-1,0); - v3d[i-1][1] = vtxCoords(elm2VtxConn(elm,i)-1,1); - v3d[i-1][2] = vtxCoords(elm2VtxConn(elm,i)-1,2); - } - v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); - v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); - v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2); - - double basisByArea3d[maxVtxsPerElm] = {0.0}; - initArray(basisByArea3d,maxVtxsPerElm,0.0); - - // calc basis - getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d); - - // interpolation step - for(int entry=0; entryparallel_for(interpolation, "interpolation"); - pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds()); -} - - -inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){ - Kokkos::Timer timer; - auto p_mesh = mpMesh.p_mesh; - auto vtxCoords = p_mesh->getMeshField(); - int numVtxs = p_mesh->getNumVertices(); - auto elm2VtxConn = p_mesh->getElm2VtxConn(); - - auto p_MPs = mpMesh.p_MPs; - auto MPsPosition = p_MPs->getPositions(); - double radius = p_mesh->getSphereRadius(); - PMT_ALWAYS_ASSERT(radius > 0); - - constexpr MeshFieldIndex meshFieldIndex1 = polyMPO::MeshF_RotLatLonIncr; - constexpr MeshFieldIndex meshFieldIndex2 = polyMPO::MeshF_OnSurfVeloIncr; - - auto meshField1 = p_mesh->getMeshField(); - auto meshField2 = p_mesh->getMeshField(); - - constexpr MaterialPointSlice mpfIndex1 = meshFieldIndexToMPSlice; - constexpr MaterialPointSlice mpfIndex2 = meshFieldIndexToMPSlice; - - const int numEntries1 = mpSliceToNumEntries(); - const int numEntries2 = mpSliceToNumEntries(); - - auto mpField1 = p_MPs->getData(); - auto mpField2 = p_MPs->getData(); - - auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) { - if(mask) { - Vec3d position3d(MPsPosition(mp, 0), MPsPosition(mp, 1), MPsPosition(mp, 2)); - Vec3d v3d[maxVtxsPerElm + 1]; - int numVtx = elm2VtxConn(elm, 0); - for (int i = 1; i <= numVtx; i++) { - v3d[i-1][0] = vtxCoords(elm2VtxConn(elm, i) - 1, 0); - v3d[i-1][1] = vtxCoords(elm2VtxConn(elm, i) - 1, 1); - v3d[i-1][2] = vtxCoords(elm2VtxConn(elm, i) - 1, 2); - } - v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0); - v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1); - v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2); - - double basisByArea3d[maxVtxsPerElm] = {0.0}; - initArray(basisByArea3d, maxVtxsPerElm, 0.0); - - getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d); - - for(int entry=0; entryparallel_for(interpolation, "sphericalInterpolationMultiField"); - pumipic::RecordTime("PolyMPO_sphericalInterpolationDispVelIncr", timer.seconds()); - } - - } //namespace polyMPO end #endif diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 292d806..4ac1c6c 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -14,9 +14,9 @@ subroutine setProcWedges(mpMesh, nCells, comm_size, nEdgesOnCell, verticesOnCell integer, dimension(:,:), pointer :: verticesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: lonCell real(kind=MPAS_RKIND) :: normalizedLong, min, max - + allocate(owningProc(nCells)) - + min = 1000 max = -1000 do i = 1, nCells @@ -74,10 +74,8 @@ subroutine runAdvectionTest2(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius, -numPush) call polympo_push(mpMesh) - - end subroutine - + end subroutine subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2Elm, & latVertex, lonVertex, nEdgesOnCell, verticesOnCell, sphereRadius) @@ -108,14 +106,17 @@ subroutine runReconstructionTest(mpMesh, numMPs, numPush, nCells, nVertices, mp2 call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass)) call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel)) - - ! Test push reconstruction + ! Although this test just does 0th order reconstruction testing, and just needs the BasisSlice, + ! calculating the coefficeints too as that will involve calculating the Basis Slice + call polympo_reconstruct_coeff_with_MPI(mpmesh) + ! Test push reconstruction do j = 1, numPush call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFElmType()) call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFVtxType()) call polympo_setReconstructionOfVel(mpMesh,0,polympo_getMeshFVtxType()) + call polympo_applyReconstruction(mpMesh) call polympo_push(mpMesh) call polympo_getMeshElmMass(mpMesh,nCells,c_loc(meshElmMass)) call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass)) @@ -163,7 +164,7 @@ subroutine runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPo meshVtxMass = TEST_VAL meshElmMass = TEST_VAL meshVtxVel = TEST_VAL - + do j = 1, numPush call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition)) call polympo_setMeshVtxOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr)) @@ -257,6 +258,7 @@ program main verticesOnCell, cellsOnCell) if (onSphere .ne. 'YES') then write (*,*) "The mesh is not spherical!" + write (*,*) "Gnomonic Projection is used currently for spehrical meshes, alternative for planar meshes not supported !" call exit(1) end if @@ -271,6 +273,8 @@ program main xCell, yCell, zCell, & verticesOnCell, cellsOnCell) + call polympo_setGnomonicProjection(mpMesh) + call polympo_setMPICommunicator(mpMesh, mpi_comm_handle); !createMPs @@ -362,7 +366,7 @@ program main call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive)) call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon)) call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition)) - + if (testType == "API") then call runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPosition, xVertex, yVertex, zVertex, latVertex) else if (testType == "MIGRATION") then diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90 index 4a9e910..c98f1b8 100644 --- a/test/testFortranMPReconstruction.f90 +++ b/test/testFortranMPReconstruction.f90 @@ -29,7 +29,7 @@ program main real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell real(kind=MPAS_RKIND), dimension(:), pointer :: areaTriangle - + integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell integer :: numMPs, vID integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, globalElms @@ -64,6 +64,7 @@ program main verticesOnCell, cellsOnCell, areaTriangle) if (onSphere .ne. 'YES') then write (*,*) "The mesh is not spherical!" + write (*,*) "Gnomonic Projection is used currently for spehrical meshes, alternative for planar meshes not supported !" call exit(1) end if @@ -78,7 +79,8 @@ program main latVertex, & xCell, yCell, zCell, & verticesOnCell, cellsOnCell, areaTriangle) - + + call polympo_setGnomonicProjection(mpMesh) nCompsDisp = 2 allocate(dispIncr(nCompsDisp,nVertices)) !createMPs @@ -110,7 +112,7 @@ program main mpPosition(2,i) = yCell(i) mpPosition(3,i) = zCell(i) end do - + call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms)) call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive)) call polympo_setMPICommunicator(mpMesh, mpi_comm_handle) @@ -120,17 +122,10 @@ program main call polympo_setMPMass(mpMesh,1,numMPs,c_loc(mpMass)) call polympo_setMPVel(mpMesh,2,numMPs,c_loc(mpVel)) - ! Test vtx reconstruction - call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFVtxType()) - call polympo_applyReconstruction(mpMesh) - call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass)) - - do i = 1, nVertices - call assert(meshVtxMass(i) < TEST_VAL+TOLERANCE .and. meshVtxMass(i) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass") - !write(*, *) 'The value of LRV is:', meshVtxMass(i) - end do - + !First Order reconstruction done before 0th order reconstruction as calculation of coefficients will + !fill the MpBasis slice !Test vtx order 1 reconstruction + call polympo_reconstruct_coeff_with_MPI(mpmesh) call polympo_setReconstructionOfMass(mpMesh,1,polympo_getMeshFVtxType()) call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType()) call polympo_applyReconstruction(mpMesh) @@ -138,9 +133,17 @@ program main call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv)) do i = 1, nVertices call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass order 1") - call assert(meshVtxVelu(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1") - call assert(meshVtxVelv(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1") - write(*, *) 'The value of LRV is:', meshVtxVelu(i)-TEST_VAL, meshVtxVelv(i)-TEST_VAL + call assert(meshVtxVelu(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1") + call assert(meshVtxVelv(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1") + end do + + ! Test vtx order 0 reconstruction + call polympo_setReconstructionOfMass(mpMesh,0,polympo_getMeshFVtxType()) + call polympo_applyReconstruction(mpMesh) + call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass)) + + do i = 1, nVertices + call assert(meshVtxMass(i) < TEST_VAL+TOLERANCE .and. meshVtxMass(i) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass") end do @@ -162,8 +165,7 @@ program main .and. meshElmMass(mp2Elm(i)) > TEST_VAL-TOLERANCE, "Error: wrong elm mass") end do end do - - + call polympo_deleteMPMesh(mpMesh) call polympo_finalize()