From d1e4e7f64b07b2ae2516a453ed961fac2749c2a6 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Thu, 8 Feb 2024 19:57:41 -0500 Subject: [PATCH 01/26] migration test --- test/CMakeLists.txt | 2 ++ test/testMigration.cpp | 22 ++++++++++++++++++++++ 2 files changed, 24 insertions(+) create mode 100644 test/testMigration.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 637a5b9f..250eb3ad 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -48,6 +48,7 @@ pmpo_add_exe(testReconstruction testReconstruction.cpp) pmpo_add_exe(testTracking testTracking.cpp) pmpo_add_exe(testRebuild testRebuild.cpp) pmpo_add_exe(timeAssmblyWachspress testTiming.cpp) +pmpo_add_exe(testMigration testMigration.cpp) pmpo_add_fortran_exe(testFortranInit testFortranInit.f90) pmpo_add_fortran_exe(testFortran testFortran.f90) pmpo_add_fortran_exe(testFortranReadMPAS testFortranReadMPAS.f90) @@ -95,3 +96,4 @@ pmpo_add_test(testFortranCreateRebuildMPs ./testFortranCreateRebuildMPs ${TEST_N pmpo_add_test(test_timing ./timeAssmblyWachspress 1 ${TEST_NC_FILE}) pmpo_add_test(test_wachspress ./testWachspress ${TEST_NC_FILE}) +pmpo_add_test(test_migration ./testMigration ${TEST_NC_FILE}) \ No newline at end of file diff --git a/test/testMigration.cpp b/test/testMigration.cpp new file mode 100644 index 00000000..df20696e --- /dev/null +++ b/test/testMigration.cpp @@ -0,0 +1,22 @@ +#include "pmpo_MPMesh.hpp" +#include "testUtils.hpp" + +using namespace polyMPO; + +int main(int argc, char* argv[] ) { + MPI_Init(&argc, &argv); + Kokkos::initialize(argc,argv); + + { + Mesh* mesh = NULL; + MPMesh* mpmesh = NULL; + void* p; + setWithMPASMeshByFortran(&p, argv[1], (int)strlen(argv[1])); + mpmesh = (MPMesh*)p; + mesh = mpmesh->p_mesh; + } + + Kokkos::finalize(); + MPI_Finalize(); + return 0; +} \ No newline at end of file From 41eef0ad84256685ae4bf8a920b34537f5862723 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Thu, 8 Feb 2024 20:19:25 -0500 Subject: [PATCH 02/26] copied old code --- src/pmpo_MPMesh.cpp | 35 ++++++++++++++++++++++++++++------- src/pmpo_materialPoints.cpp | 19 +++++++++++++++++++ src/pmpo_materialPoints.hpp | 14 ++++++++++---- src/pmpo_mesh.cpp | 5 +++++ src/pmpo_mesh.hpp | 2 ++ 5 files changed, 64 insertions(+), 11 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 322317f8..8a4e451b 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -84,7 +84,9 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ auto mpPositions = p_MPs->getData(); auto mpTgtPos = p_MPs->getData(); - auto MPs2Elm = p_MPs->getData();; + auto MPs2Elm = p_MPs->getData(); + auto MPs2Proc = p_MPs->getData(); + auto elm2Process = p_mesh->getElm2Process(); if(printVTPIndex>=0) { printVTP_mesh(printVTPIndex); @@ -132,6 +134,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ } if(closestElm<0){ MPs2Elm(mp) = iElm; + MPs2Proc(mp) = elm2Process(iElm); break; }else{ iElm = closestElm; @@ -258,17 +261,35 @@ void MPMesh::T2LTracking(Vec2dView dx){ p_MPs->parallel_for(T2LCalc,"T2lTrackingCalc"); } +bool getAnyIsMigrating(bool isMigrating) { + int comm_rank; + MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank); + int comm_size; + MPI_Comm_size(MPI_COMM_WORLD, &comm_size); + + bool anyIsMigrating = isMigrating; + for (int i=0; i < comm_size; i++) { + if (i == comm_rank) MPI_Bcast(&anyIsMigrating, 1, MPI_C_BOOL, i, MPI_COMM_WORLD); + else MPI_Bcast(&isMigrating, 1, MPI_C_BOOL, i, MPI_COMM_WORLD); + anyIsMigrating |= isMigrating; + } + return anyIsMigrating; +} + void MPMesh::push(){ p_mesh->computeRotLatLonIncr(); sphericalInterpolation(*this); p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ - CVTTrackingElmCenterBased(); // move to Tgt_XYZ - - p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ - p_MPs->updateMPSlice(); // Tgt becomes Cur - p_MPs->rebuild(); //rebuild pumi-pic - p_MPs->updateMPElmID(); //update mpElm IDs slices + while(true) { + CVTTrackingElmCenterBased(); // move to Tgt_XYZ + p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ + p_MPs->updateMPSlice(); // Tgt becomes Cur + bool isMigrating = p_MPs->migrate(); + bool anyIsMigrating = getAnyIsMigrating(isMigrating); + p_MPs->updateMPElmID(); //update mpElm IDs slices + if (!anyIsMigrating) break; + } } void MPMesh::printVTP_mesh(int printVTPIndex){ diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index c0c53b69..08219a04 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -101,6 +101,25 @@ void MaterialPoints::finishRebuild() { rebuildFields.ongoing = false; } +bool MaterialPoints::migrate() { + auto MPs2Elm = getData(); + auto MPs2Proc = getData(); + + IntView new_elem("new_elem", MPs->capacity()); + IntView new_process("new_process", MPs->capacity()); + IntView isMigrating("isMigrating", 1); + + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + Kokkos::parallel_for("setMigrateFields", MPs->capacity(), KOKKOS_LAMBDA(int i) { + new_elem(i) = MPs2Elm(i); + new_process(i) = MPs2Proc(i); + if (new_process(i) != rank) isMigrating(0) = 1; + }); + MPs->migrate(new_elem, new_process); + return pumipic::getLastValue(isMigrating) > 0; +} + bool MaterialPoints::rebuildOngoing() { return rebuildFields.ongoing; } void MaterialPoints::setAppIDFunc(IntFunc getAppIDIn) { getAppID = getAppIDIn; } diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp index 33beff2f..29e56ee1 100644 --- a/src/pmpo_materialPoints.hpp +++ b/src/pmpo_materialPoints.hpp @@ -21,6 +21,7 @@ using defaultSpace = Kokkos::DefaultExecutionSpace::memory_space; typedef int mp_flag_t; typedef int mp_id_t; typedef int mp_elm_id_t; +typedef int mp_proc_id_t; typedef double mp_sclr_t[1];//TODO typedef vec2d_t mp_vec2d_t; typedef double mp_vec3d_t[3];//TODO @@ -49,7 +50,8 @@ enum MaterialPointSlice { MPF_Stress_Div, MPF_Shear_Traction, MPF_Constv_Mdl_Param, - MPF_MP_APP_ID + MPF_MP_APP_ID, + MPF_Tgt_Proc_ID }; enum Operating_Mode{ @@ -76,7 +78,8 @@ const static std::map> {MPF_Stress_Div, {3,MeshF_Unsupported}}, {MPF_Shear_Traction, {3,MeshF_Unsupported}}, {MPF_Constv_Mdl_Param,{12,MeshF_Unsupported}}, - {MPF_MP_APP_ID, {1,MeshF_Invalid}}}; + {MPF_MP_APP_ID, {1,MeshF_Invalid}}, + {MPF_Tgt_Proc_ID, {0,MeshF_Invalid}}}; //TODO: What does this integer mean? const static std::vector> mpSliceSwap = {{MPF_Cur_Elm_ID, MPF_Tgt_Elm_ID}, @@ -101,7 +104,8 @@ typedef MemberTypesMaterialPointTypes; typedef ps::ParticleStructure PS; @@ -134,7 +138,9 @@ class MaterialPoints { void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID, Kokkos::View addedMPMask); void finishRebuild(); bool rebuildOngoing(); - + + bool migrate(); + template typename std::enable_if::type setRebuildMPSlice(mpSliceData mpSliceIn); diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index 78b3a7be..70c4d5f4 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -55,4 +55,9 @@ namespace polyMPO{ }); } + IntView Mesh::getElm2Process() { + IntView elm2Process("elm2Process", numElms_); + return elm2Process; + } + } // namespace polyMPO diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 3f94256a..9e6aa92c 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -103,6 +103,8 @@ class Mesh { bool checkMeshType(int meshType); bool checkGeomType(int geomType); + IntView getElm2Process(); + mesh_type getMeshType() { return meshType_; } geom_type getGeomType() { return geomType_; } double getSphereRadius() { return sphereRadius_; } From 8a79fa44699c25d9b3c745845ee60a5fac1b2b65 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Fri, 9 Feb 2024 15:35:02 -0500 Subject: [PATCH 03/26] owning proc --- src/pmpo_c.cpp | 14 ++++++++++++++ src/pmpo_c.h | 1 + src/pmpo_fortran.f90 | 13 +++++++++++++ src/pmpo_mesh.hpp | 7 ++++++- test/testMigration.cpp | 25 +++++++++++++++++++++++++ 5 files changed, 59 insertions(+), 1 deletion(-) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index df280d48..c447dd6d 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -642,3 +642,17 @@ void polympo_push_f(MPMesh_ptr p_mpmesh){ checkMPMeshValid(p_mpmesh); ((polyMPO::MPMesh*)p_mpmesh) ->push(); } + +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); +} \ No newline at end of file diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 23200b42..dfd6116b 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -50,6 +50,7 @@ void polympo_getMeshNumElms_f(MPMesh_ptr p_mpmesh, int & numElms); void polympo_setMeshNumEdgesPerElm_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array); 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); //Mesh fields void polympo_setMeshVtxCoords_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* xArray, const double* yArray, const double* zArray); diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index d2051ee3..3accf1c1 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -507,6 +507,19 @@ subroutine polympo_getMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) & type(c_ptr), value :: array end subroutine !--------------------------------------------------------------------------- + !> @brief set the index to cell id array from a host 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 + !--------------------------------------------------------------------------- !> @brief calculate the MPs from given mesh vertices rotational latitude !> longitude, update the MP slices !> MPs MUST have rotated flag set to True(>0) diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 9e6aa92c..1ccfd5cc 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -66,7 +66,9 @@ class Mesh { //IntView nEdgesPerElm_; IntVtx2ElmView elm2VtxConn_; IntElm2ElmView elm2ElmConn_; - + + IntView owningProc_; + //start of meshFields DoubleVec3dView vtxCoords_; DoubleSclrView vtxRotLat_; @@ -112,6 +114,7 @@ class Mesh { int getNumElements() { return numElms_; } IntVtx2ElmView getElm2VtxConn() { return elm2VtxConn_; } IntElm2ElmView getElm2ElmConn() { return elm2ElmConn_; } + IntView getOnwningProc() { return owningProc_; } template auto getMeshField(); void setMeshVtxBasedFieldSize(); @@ -131,6 +134,8 @@ 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 computeRotLatLonIncr(); }; diff --git a/test/testMigration.cpp b/test/testMigration.cpp index df20696e..f6e3fb76 100644 --- a/test/testMigration.cpp +++ b/test/testMigration.cpp @@ -14,6 +14,31 @@ int main(int argc, char* argv[] ) { setWithMPASMeshByFortran(&p, argv[1], (int)strlen(argv[1])); mpmesh = (MPMesh*)p; mesh = mpmesh->p_mesh; + + int numElms = mesh->getNumElements(); + auto elm2VtxConn = mesh->getElm2VtxConn(); + auto vtxCoords = mesh->getMeshField(); + + Vec3dView elmCenter("elementCenter",numElms); + Kokkos::parallel_for("calcElementCenter", numElms, KOKKOS_LAMBDA(const int elm){ + int numVtx = elm2VtxConn(elm,0); + double sum_x = 0.0, sum_y = 0.0, sum_z = 0.0; + for(int i=1; i<= numVtx; i++){ + sum_x += vtxCoords(elm2VtxConn(elm,i)-1,0); + sum_y += vtxCoords(elm2VtxConn(elm,i)-1,1); + sum_z += vtxCoords(elm2VtxConn(elm,i)-1,2); + } + elmCenter(elm)[0] = sum_x/numVtx; + elmCenter(elm)[1] = sum_y/numVtx; + elmCenter(elm)[2] = sum_z/numVtx; + }); + + IntView owningProc("owningProc", numElms); + Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){ + owningProc(elm) = elmCenter(elm)[1] > 0 ? 1 : 0; + }); + mesh->setMeshEdit(true); + mesh->setOwningProc(owningProc); } Kokkos::finalize(); From 8234bc2fd8c3520f843befc51d4a53c8df9932f1 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Fri, 9 Feb 2024 17:09:03 -0500 Subject: [PATCH 04/26] initialize material points --- test/testMigration.cpp | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/test/testMigration.cpp b/test/testMigration.cpp index f6e3fb76..5eaf5111 100644 --- a/test/testMigration.cpp +++ b/test/testMigration.cpp @@ -8,14 +8,28 @@ int main(int argc, char* argv[] ) { Kokkos::initialize(argc,argv); { - Mesh* mesh = NULL; MPMesh* mpmesh = NULL; void* p; setWithMPASMeshByFortran(&p, argv[1], (int)strlen(argv[1])); mpmesh = (MPMesh*)p; - mesh = mpmesh->p_mesh; - + Mesh* mesh = mpmesh->p_mesh; int numElms = mesh->getNumElements(); + int numMPs = numElms; + + IntView mpsPerElm("mpsPerElm", numElms); + IntView activeMP2Elm("activeMP2Elm", numMPs); + IntView activeMPIDs("activeMPIDs", numMPs); + + Kokkos::parallel_for("setMPsPerElm", numElms, KOKKOS_LAMBDA(const int elm){ + mpsPerElm(elm) = 1; + }); + Kokkos::parallel_for("setMPs", numMPs, KOKKOS_LAMBDA(const int mp){ + activeMP2Elm(mp) = 1; + activeMPIDs(mp) = mp; + }); + + MaterialPoints p_MPs(numElms, numMPs, mpsPerElm, activeMP2Elm, activeMPIDs); + auto elm2VtxConn = mesh->getElm2VtxConn(); auto vtxCoords = mesh->getMeshField(); From 4bad5b8746db4c5709f6e38279eb887d967d44d6 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Tue, 13 Feb 2024 16:00:11 -0500 Subject: [PATCH 05/26] migration test working --- src/pmpo_MPMesh.cpp | 8 +++++--- src/pmpo_materialPoints.cpp | 13 ++++++++----- src/pmpo_mesh.cpp | 3 +-- src/pmpo_mesh.hpp | 1 - test/CMakeLists.txt | 34 +++++++++++++++++----------------- test/testMigration.cpp | 20 ++++++++++++++++++-- 6 files changed, 49 insertions(+), 30 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 8a4e451b..70b05f62 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -134,7 +134,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ } if(closestElm<0){ MPs2Elm(mp) = iElm; - MPs2Proc(mp) = elm2Process(iElm); + if (elm2Process.size() > 0) + MPs2Proc(mp) = elm2Process(iElm); break; }else{ iElm = closestElm; @@ -283,13 +284,14 @@ void MPMesh::push(){ while(true) { CVTTrackingElmCenterBased(); // move to Tgt_XYZ - p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ - p_MPs->updateMPSlice(); // Tgt becomes Cur bool isMigrating = p_MPs->migrate(); bool anyIsMigrating = getAnyIsMigrating(isMigrating); p_MPs->updateMPElmID(); //update mpElm IDs slices if (!anyIsMigrating) break; } + + p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ + p_MPs->updateMPSlice(); // Tgt becomes Cur } void MPMesh::printVTP_mesh(int printVTPIndex){ diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index 08219a04..b044256c 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -111,11 +111,14 @@ bool MaterialPoints::migrate() { int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); - Kokkos::parallel_for("setMigrateFields", MPs->capacity(), KOKKOS_LAMBDA(int i) { - new_elem(i) = MPs2Elm(i); - new_process(i) = MPs2Proc(i); - if (new_process(i) != rank) isMigrating(0) = 1; - }); + auto setMigrationFields = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { + if (mask) { + new_elem(mp) = MPs2Elm(mp); + new_process(mp) = MPs2Proc(mp); + if (new_process(mp) != rank) isMigrating(0) = 1; + } + }; + parallel_for(setMigrationFields, "setMigrationFields"); MPs->migrate(new_elem, new_process); return pumipic::getLastValue(isMigrating) > 0; } diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index 70c4d5f4..084e4352 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -56,8 +56,7 @@ namespace polyMPO{ } IntView Mesh::getElm2Process() { - IntView elm2Process("elm2Process", numElms_); - return elm2Process; + return owningProc_; } } // namespace polyMPO diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 1ccfd5cc..22e4d8b2 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -114,7 +114,6 @@ class Mesh { int getNumElements() { return numElms_; } IntVtx2ElmView getElm2VtxConn() { return elm2VtxConn_; } IntElm2ElmView getElm2ElmConn() { return elm2ElmConn_; } - IntView getOnwningProc() { return owningProc_; } template auto getMeshField(); void setMeshVtxBasedFieldSize(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 250eb3ad..3eaa813a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -60,24 +60,24 @@ pmpo_add_fortran_exe(testFortranCreateRebuildMPs testFortranCreateRebuildMPs.f90 #add test find_program(MPIRUN_EXECUTABLE NAMES mpirun) -function(pmpo_add_test testname command) +function(pmpo_add_test testname procs command) if(MPIRUN_EXECUTABLE) add_test(NAME ${testname} - COMMAND mpirun --bind-to core -np 1 ${command} ${ARGN}) + COMMAND mpirun --bind-to core -np ${procs} ${command} ${ARGN}) else() add_test(NAME ${testname} COMMAND ${command} ${ARGN}) endif() endfunction() -pmpo_add_test(unit_test ./unitTest) -pmpo_add_test(mp_unit_test ./mpUnitTest) -pmpo_add_test(test_reconstruction ./testReconstruction) -pmpo_add_test(test_tracking ./testTracking) -pmpo_add_test(test_rebuild ./testRebuild) -pmpo_add_test(testFortranInit ./testFortranInit) -pmpo_add_test(testFortran ./testFortran) -pmpo_add_test(testFortranMPMeshModule ./testFortranMPMeshModule) -pmpo_add_test(testFortranMPAppIDs ./testFortranMPAppIDs) +pmpo_add_test(unit_test 1 ./unitTest) +pmpo_add_test(mp_unit_test 1 ./mpUnitTest) +pmpo_add_test(test_reconstruction 1 ./testReconstruction) +pmpo_add_test(test_tracking 1 ./testTracking) +pmpo_add_test(test_rebuild 1 ./testRebuild) +pmpo_add_test(testFortranInit 1 ./testFortranInit) +pmpo_add_test(testFortran 1 ./testFortran) +pmpo_add_test(testFortranMPMeshModule 1 ./testFortranMPMeshModule) +pmpo_add_test(testFortranMPAppIDs 1 ./testFortranMPAppIDs) #set NC file for test set(TEST_NC_FILE_PLANAR "${CMAKE_SOURCE_DIR}/test/sample_mpas_meshes/planar_nonuniform_cvt_for_square_673elms.nc") @@ -87,13 +87,13 @@ if(EXISTS ${TEST_NC_FILE_SPHERICAL}) message(STATUS "Found Spherical NC File: ${TEST_NC_FILE_SPHERICAL}") set(TEST_NC_FILE ${TEST_NC_FILE_SPHERICAL}) # this test only works for spherical NC files! - pmpo_add_test(testFortranMPAdvection ./testFortranMPAdvection ${TEST_NC_FILE}) + pmpo_add_test(testFortranMPAdvection 1 ./testFortranMPAdvection ${TEST_NC_FILE}) else() set(TEST_NC_FILE ${TEST_NC_FILE_PLANAR}) endif() -pmpo_add_test(testFortranReadMPAS ./testFortranReadMPAS ${TEST_NC_FILE}) -pmpo_add_test(testFortranCreateRebuildMPs ./testFortranCreateRebuildMPs ${TEST_NC_FILE}) +pmpo_add_test(testFortranReadMPAS 1 ./testFortranReadMPAS ${TEST_NC_FILE}) +pmpo_add_test(testFortranCreateRebuildMPs 1 ./testFortranCreateRebuildMPs ${TEST_NC_FILE}) -pmpo_add_test(test_timing ./timeAssmblyWachspress 1 ${TEST_NC_FILE}) -pmpo_add_test(test_wachspress ./testWachspress ${TEST_NC_FILE}) -pmpo_add_test(test_migration ./testMigration ${TEST_NC_FILE}) \ No newline at end of file +pmpo_add_test(test_timing 1 ./timeAssmblyWachspress 1 ${TEST_NC_FILE}) +pmpo_add_test(test_wachspress 1 ./testWachspress ${TEST_NC_FILE}) +pmpo_add_test(test_migration 2 ./testMigration ${TEST_NC_FILE}) \ No newline at end of file diff --git a/test/testMigration.cpp b/test/testMigration.cpp index 5eaf5111..9995c941 100644 --- a/test/testMigration.cpp +++ b/test/testMigration.cpp @@ -7,6 +7,11 @@ int main(int argc, char* argv[] ) { MPI_Init(&argc, &argv); Kokkos::initialize(argc,argv); + int comm_rank; + MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank); + int comm_size; + MPI_Comm_size(MPI_COMM_WORLD, &comm_size); + { MPMesh* mpmesh = NULL; void* p; @@ -24,11 +29,12 @@ int main(int argc, char* argv[] ) { mpsPerElm(elm) = 1; }); Kokkos::parallel_for("setMPs", numMPs, KOKKOS_LAMBDA(const int mp){ - activeMP2Elm(mp) = 1; + activeMP2Elm(mp) = mp; activeMPIDs(mp) = mp; }); MaterialPoints p_MPs(numElms, numMPs, mpsPerElm, activeMP2Elm, activeMPIDs); + mpmesh->p_MPs = &p_MPs; auto elm2VtxConn = mesh->getElm2VtxConn(); auto vtxCoords = mesh->getMeshField(); @@ -49,10 +55,20 @@ int main(int argc, char* argv[] ) { IntView owningProc("owningProc", numElms); Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){ - owningProc(elm) = elmCenter(elm)[1] > 0 ? 1 : 0; + owningProc(elm) = elmCenter(elm)[1] > 0 ? comm_size-1 : 0; }); mesh->setMeshEdit(true); mesh->setOwningProc(owningProc); + mpmesh->push(); + + auto MPs2Proc = p_MPs.getData(); + auto checkPostBackMigrate = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { + if (mask) { + assert(MPs2Proc(mp) == comm_rank); + assert(owningProc(e) == comm_rank); + } + }; + p_MPs.parallel_for(checkPostBackMigrate, "checkPostBackMigrate"); } Kokkos::finalize(); From 007c14e761b76d7ee7ef25e1a931900b4120f98e Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 14 Feb 2024 15:21:55 -0500 Subject: [PATCH 06/26] moved out divide procs --- test/testMigration.cpp | 47 +++++++++++++++++++++++------------------- 1 file changed, 26 insertions(+), 21 deletions(-) diff --git a/test/testMigration.cpp b/test/testMigration.cpp index 9995c941..4b026ca6 100644 --- a/test/testMigration.cpp +++ b/test/testMigration.cpp @@ -3,6 +3,31 @@ using namespace polyMPO; +IntView divideProcsInHalf(Mesh* mesh, int numElms, int comm_rank, int comm_size) { + auto elm2VtxConn = mesh->getElm2VtxConn(); + auto vtxCoords = mesh->getMeshField(); + + Vec3dView elmCenter("elementCenter",numElms); + Kokkos::parallel_for("calcElementCenter", numElms, KOKKOS_LAMBDA(const int elm){ + int numVtx = elm2VtxConn(elm,0); + double sum_x = 0.0, sum_y = 0.0, sum_z = 0.0; + for(int i=1; i<= numVtx; i++){ + sum_x += vtxCoords(elm2VtxConn(elm,i)-1,0); + sum_y += vtxCoords(elm2VtxConn(elm,i)-1,1); + sum_z += vtxCoords(elm2VtxConn(elm,i)-1,2); + } + elmCenter(elm)[0] = sum_x/numVtx; + elmCenter(elm)[1] = sum_y/numVtx; + elmCenter(elm)[2] = sum_z/numVtx; + }); + + IntView owningProc("owningProc", numElms); + Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){ + owningProc(elm) = elmCenter(elm)[1] > 0 ? comm_size-1 : 0; + }); + return owningProc; +} + int main(int argc, char* argv[] ) { MPI_Init(&argc, &argv); Kokkos::initialize(argc,argv); @@ -36,27 +61,7 @@ int main(int argc, char* argv[] ) { MaterialPoints p_MPs(numElms, numMPs, mpsPerElm, activeMP2Elm, activeMPIDs); mpmesh->p_MPs = &p_MPs; - auto elm2VtxConn = mesh->getElm2VtxConn(); - auto vtxCoords = mesh->getMeshField(); - - Vec3dView elmCenter("elementCenter",numElms); - Kokkos::parallel_for("calcElementCenter", numElms, KOKKOS_LAMBDA(const int elm){ - int numVtx = elm2VtxConn(elm,0); - double sum_x = 0.0, sum_y = 0.0, sum_z = 0.0; - for(int i=1; i<= numVtx; i++){ - sum_x += vtxCoords(elm2VtxConn(elm,i)-1,0); - sum_y += vtxCoords(elm2VtxConn(elm,i)-1,1); - sum_z += vtxCoords(elm2VtxConn(elm,i)-1,2); - } - elmCenter(elm)[0] = sum_x/numVtx; - elmCenter(elm)[1] = sum_y/numVtx; - elmCenter(elm)[2] = sum_z/numVtx; - }); - - IntView owningProc("owningProc", numElms); - Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){ - owningProc(elm) = elmCenter(elm)[1] > 0 ? comm_size-1 : 0; - }); + IntView owningProc = divideProcsInHalf(mesh, numElms, comm_rank, comm_size); mesh->setMeshEdit(true); mesh->setOwningProc(owningProc); mpmesh->push(); From 851f99c678b4daf213bcb8189f8157d8e4dd1edf Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 14 Feb 2024 17:33:20 -0500 Subject: [PATCH 07/26] divide in wedges --- test/CMakeLists.txt | 2 +- test/testMigration.cpp | 27 ++++++++++++++------------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3eaa813a..4fe5cd5f 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -96,4 +96,4 @@ pmpo_add_test(testFortranCreateRebuildMPs 1 ./testFortranCreateRebuildMPs ${TEST pmpo_add_test(test_timing 1 ./timeAssmblyWachspress 1 ${TEST_NC_FILE}) pmpo_add_test(test_wachspress 1 ./testWachspress ${TEST_NC_FILE}) -pmpo_add_test(test_migration 2 ./testMigration ${TEST_NC_FILE}) \ No newline at end of file +pmpo_add_test(test_migration 8 ./testMigration ${TEST_NC_FILE}) \ No newline at end of file diff --git a/test/testMigration.cpp b/test/testMigration.cpp index 4b026ca6..11f993c7 100644 --- a/test/testMigration.cpp +++ b/test/testMigration.cpp @@ -3,27 +3,28 @@ using namespace polyMPO; -IntView divideProcsInHalf(Mesh* mesh, int numElms, int comm_rank, int comm_size) { +IntView divideProcsInWedges(Mesh* mesh, int numElms, int comm_rank, int comm_size) { auto elm2VtxConn = mesh->getElm2VtxConn(); - auto vtxCoords = mesh->getMeshField(); + auto vtxLat = mesh->getMeshField(); - Vec3dView elmCenter("elementCenter",numElms); - Kokkos::parallel_for("calcElementCenter", numElms, KOKKOS_LAMBDA(const int elm){ + DoubleView elmLat("elmLat",numElms); + DoubleView max("max", 1); + DoubleView min("min", 1); + Kokkos::parallel_for("calcElementLat", numElms, KOKKOS_LAMBDA(const int elm){ int numVtx = elm2VtxConn(elm,0); - double sum_x = 0.0, sum_y = 0.0, sum_z = 0.0; + double sum = 0.0; for(int i=1; i<= numVtx; i++){ - sum_x += vtxCoords(elm2VtxConn(elm,i)-1,0); - sum_y += vtxCoords(elm2VtxConn(elm,i)-1,1); - sum_z += vtxCoords(elm2VtxConn(elm,i)-1,2); + sum += vtxLat(elm2VtxConn(elm,i)-1); } - elmCenter(elm)[0] = sum_x/numVtx; - elmCenter(elm)[1] = sum_y/numVtx; - elmCenter(elm)[2] = sum_z/numVtx; + elmLat(elm) = sum/numVtx; + Kokkos::atomic_max(&max(0), elmLat(elm)); + Kokkos::atomic_min(&min(0), elmLat(elm)); }); IntView owningProc("owningProc", numElms); Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){ - owningProc(elm) = elmCenter(elm)[1] > 0 ? comm_size-1 : 0; + double normalizedLat = (elmLat(elm) - min(0)) / (max(0) - min(0)); + owningProc(elm) = normalizedLat * comm_size; }); return owningProc; } @@ -61,7 +62,7 @@ int main(int argc, char* argv[] ) { MaterialPoints p_MPs(numElms, numMPs, mpsPerElm, activeMP2Elm, activeMPIDs); mpmesh->p_MPs = &p_MPs; - IntView owningProc = divideProcsInHalf(mesh, numElms, comm_rank, comm_size); + IntView owningProc = divideProcsInWedges(mesh, numElms, comm_rank, comm_size); mesh->setMeshEdit(true); mesh->setOwningProc(owningProc); mpmesh->push(); From 66e00fddaba822c2f683972559011773eafea524 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 14 Feb 2024 17:48:58 -0500 Subject: [PATCH 08/26] rename procWedges --- test/testMigration.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/testMigration.cpp b/test/testMigration.cpp index 11f993c7..bba49067 100644 --- a/test/testMigration.cpp +++ b/test/testMigration.cpp @@ -3,7 +3,7 @@ using namespace polyMPO; -IntView divideProcsInWedges(Mesh* mesh, int numElms, int comm_rank, int comm_size) { +IntView procWedges(Mesh* mesh, int numElms, int comm_size) { auto elm2VtxConn = mesh->getElm2VtxConn(); auto vtxLat = mesh->getMeshField(); @@ -62,7 +62,7 @@ int main(int argc, char* argv[] ) { MaterialPoints p_MPs(numElms, numMPs, mpsPerElm, activeMP2Elm, activeMPIDs); mpmesh->p_MPs = &p_MPs; - IntView owningProc = divideProcsInWedges(mesh, numElms, comm_rank, comm_size); + IntView owningProc = procWedges(mesh, numElms, comm_size); mesh->setMeshEdit(true); mesh->setOwningProc(owningProc); mpmesh->push(); From 561ee392c5c9e251b675864b74e646ef71777bda Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 21 Feb 2024 14:00:40 -0500 Subject: [PATCH 09/26] use mpi_allreduce --- src/pmpo_MPMesh.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 70b05f62..e8379342 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -268,12 +268,8 @@ bool getAnyIsMigrating(bool isMigrating) { int comm_size; MPI_Comm_size(MPI_COMM_WORLD, &comm_size); - bool anyIsMigrating = isMigrating; - for (int i=0; i < comm_size; i++) { - if (i == comm_rank) MPI_Bcast(&anyIsMigrating, 1, MPI_C_BOOL, i, MPI_COMM_WORLD); - else MPI_Bcast(&isMigrating, 1, MPI_C_BOOL, i, MPI_COMM_WORLD); - anyIsMigrating |= isMigrating; - } + bool anyIsMigrating = false; + MPI_Allreduce(&isMigrating, &anyIsMigrating, 1, MPI_C_BOOL, MPI_LOR, MPI_COMM_WORLD); return anyIsMigrating; } From 3dbf75a127738e30bc8083ad9b047f788bbb35da Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 21 Feb 2024 14:07:12 -0500 Subject: [PATCH 10/26] do while --- src/pmpo_MPMesh.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index e8379342..eed994a1 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -278,13 +278,14 @@ void MPMesh::push(){ sphericalInterpolation(*this); p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ - while(true) { + bool anyIsMigrating = false; + do { CVTTrackingElmCenterBased(); // move to Tgt_XYZ bool isMigrating = p_MPs->migrate(); - bool anyIsMigrating = getAnyIsMigrating(isMigrating); + anyIsMigrating = getAnyIsMigrating(isMigrating); p_MPs->updateMPElmID(); //update mpElm IDs slices - if (!anyIsMigrating) break; - } + } + while (anyIsMigrating); p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ p_MPs->updateMPSlice(); // Tgt becomes Cur From 59d7c6ebf00f81ff601515e95a9698759d901e1e Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 21 Feb 2024 14:08:41 -0500 Subject: [PATCH 11/26] update comment --- src/pmpo_fortran.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 3accf1c1..d7d7b9a3 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -507,7 +507,7 @@ subroutine polympo_getMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) & type(c_ptr), value :: array end subroutine !--------------------------------------------------------------------------- - !> @brief set the index to cell id array from a host array + !> @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 From 23ff539ed0a4aad852c1c1b9684cdf6289f0d50f Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 21 Feb 2024 15:07:13 -0500 Subject: [PATCH 12/26] migration timer --- src/pmpo_materialPoints.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index b044256c..a9430dc8 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -102,6 +102,7 @@ void MaterialPoints::finishRebuild() { } bool MaterialPoints::migrate() { + Kokkos::Timer timer; auto MPs2Elm = getData(); auto MPs2Proc = getData(); @@ -120,6 +121,10 @@ bool MaterialPoints::migrate() { }; parallel_for(setMigrationFields, "setMigrationFields"); MPs->migrate(new_elem, new_process); + + if (getOpMode() == polyMPO::MP_DEBUG) + printf("Material point migration: %f\n", timer.seconds()); + return pumipic::getLastValue(isMigrating) > 0; } From b7af748d182d1567a49f51ec1966b1ce65dcafef Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 21 Feb 2024 15:56:55 -0500 Subject: [PATCH 13/26] only run test if mpi --- test/CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4fe5cd5f..533a0055 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -96,4 +96,6 @@ pmpo_add_test(testFortranCreateRebuildMPs 1 ./testFortranCreateRebuildMPs ${TEST pmpo_add_test(test_timing 1 ./timeAssmblyWachspress 1 ${TEST_NC_FILE}) pmpo_add_test(test_wachspress 1 ./testWachspress ${TEST_NC_FILE}) -pmpo_add_test(test_migration 8 ./testMigration ${TEST_NC_FILE}) \ No newline at end of file +if(MPIRUN_EXECUTABLE) + pmpo_add_test(test_migration 8 ./testMigration ${TEST_NC_FILE}) +endif() \ No newline at end of file From 0ec75be413504adb31c591155ccc65e2f71c3b95 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Mon, 15 Jul 2024 20:19:23 -0700 Subject: [PATCH 14/26] added new files --- test/calculateDisplacement.f90 | 40 ++++ test/testFortranMPAdvection.f90 | 332 +++++++++++++++++--------------- 2 files changed, 213 insertions(+), 159 deletions(-) create mode 100644 test/calculateDisplacement.f90 diff --git a/test/calculateDisplacement.f90 b/test/calculateDisplacement.f90 new file mode 100644 index 00000000..dac64c2c --- /dev/null +++ b/test/calculateDisplacement.f90 @@ -0,0 +1,40 @@ +! subroutine calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) +! use :: polympo +! use :: readMPAS +! use :: iso_c_binding +! implicit none + +! real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex +! real(kind=MPAS_RKIND) :: maxlon, minlon, deltaLon, sphereRadius +! integer, dimension(:), pointer :: nEdgesOnCell +! integer, dimension(:,:), pointer :: verticesOnCell +! integer :: i, j, nVertices, nCompsDisp +! real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr +! type(c_ptr) :: mpMesh + +! nCompsDisp = 2 +! allocate(dispIncr(nCompsDisp,nVertices)) + +! ! check first element/cell for delta +! maxlon = minval(lonVertex) +! minlon = maxval(lonVertex) +! do i = 1, nEdgesOnCell(1) +! j = verticesOnCell(i,1) +! if(maxlon .lt. lonVertex(j)) then +! maxlon = lonVertex(j) +! endif +! if(minlon .gt. lonVertex(j)) then +! minlon = lonVertex(j) +! endif +! end do + +! deltaLon = maxlon - minlon + +! do i = 1,nVertices +! dispIncr(1,i) = sphereRadius*cos(latVertex(i))*deltaLon +! dispIncr(2,i) = 0.0_MPAS_RKIND +! end do +! call polympo_setMeshVtxOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr)) + +! deallocate(dispIncr) +! end subroutine \ No newline at end of file diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index a1b3fb88..6f7755fa 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -1,3 +1,29 @@ +module advectionTests + contains + include "calculateDisplacement.f90" + + subroutine runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) + use :: polympo + use :: readMPAS + use :: iso_c_binding + implicit none + + type(c_ptr) :: mpMesh + integer :: i, numPush, nVertices + real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + integer, dimension(:), pointer :: nEdgesOnCell + integer, dimension(:,:), pointer :: verticesOnCell + real(kind=MPAS_RKIND) :: sphereRadius + + + ! do i = 1, numPush + ! call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) + ! call polympo_push(mpMesh) + ! end do + + end subroutine + +end module !--------------------------------------------------------------------------- !> todo add a discription !--------------------------------------------------------------------------- @@ -5,186 +31,174 @@ program main use :: polympo use :: readMPAS use :: iso_c_binding + use :: advectionTests implicit none include 'mpif.h' !integer, parameter :: APP_RKIND = selected_real_kind(15) type(c_ptr) :: mpMesh integer :: ierr, self - integer :: argc, i, j, arglen, k + integer :: argc, i, j, arglen, k, m, mpsScaleFactorPerVtx, localNumMPs integer :: setMeshOption, setMPOption integer :: maxEdges, vertexDegree, nCells, nVertices - integer :: nCompsDisp integer :: mpi_comm_handle = MPI_COMM_WORLD - real(kind=MPAS_RKIND) :: xc, yc, zc, radius, maxlon, minlon, deltaLon, lon + real(kind=MPAS_RKIND) :: xc, yc, zc, xMP, yMP, zMP, radius, lon real(kind=MPAS_RKIND) :: pi = 4.0_MPAS_RKIND*atan(1.0_MPAS_RKIND) - character (len=2048) :: filename - real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr + character (len=2048) :: filename, input character (len=64) :: onSphere - real(kind=MPAS_RKIND) :: sphereRadius, xComputed, yComputed, zComputed, latComputed, lonComputed + real(kind=MPAS_RKIND) :: sphereRadius integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell - integer :: numMPs + integer :: numMPs, numMPsCount, numPush integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition, mpLatLon - logical :: inBound integer, parameter :: MP_ACTIVE = 1 integer, parameter :: MP_INACTIVE = 0 integer, parameter :: INVALID_ELM_ID = -1 - call mpi_init(ierr) - call mpi_comm_rank(mpi_comm_handle, self, ierr) - - call polympo_setMPICommunicator(mpi_comm_handle) - call polympo_initialize() - - call polympo_checkPrecisionForRealKind(MPAS_RKIND) - argc = command_argument_count() - if(argc == 1) then - call get_command_argument(1, filename) - else - write(0, *) "Usage: ./testFortranInterpolatePush " - end if - - call readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & - nCells, nVertices, nEdgesOnCell, & - onSphere, sphereRadius, & - xVertex, yVertex, zVertex, & - latVertex, lonVertex, & - verticesOnCell, cellsOnCell) - if (onSphere .ne. 'YES') then - write (*,*) "The mesh is not spherical!" - call exit(1) - end if - - setMeshOption = 0 !create an empty mesh - setMPOption = 0 !create an empty set of MPs - mpMesh = polympo_createMPMesh(setMeshOption,setMPOption) !creates test mesh - call loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & - nCells, nVertices, nEdgesOnCell, & - onSphere, sphereRadius, & - xVertex, yVertex, zVertex, & - latVertex, & - verticesOnCell, cellsOnCell) - - nCompsDisp = 2 - allocate(dispIncr(nCompsDisp,nVertices)) - !createMPs - numMPs = nCells - allocate(mpsPerElm(nCells)) - allocate(mp2Elm(numMPs)) - allocate(isMPActive(numMPs)) - allocate(mpPosition(3,numMPs)) - allocate(mpLatLon(2,numMPs)) - - isMPActive = MP_ACTIVE !all active MPs and some changed below - mpsPerElm = 1 !all elements have 1 MP and some changed below - do i = 1,numMPs - mp2Elm(i) = i - end do - do i = 1, numMPs - inBound = .true. - do k = 1, nEdgesOnCell(i) - j = verticesOnCell(k,i) - if ((latVertex(j) .gt. 0.4*pi) .or. (latVertex(j) .lt. -0.4*pi)) then - inBound = .false. - isMPActive(i) = MP_INACTIVE - mpsPerElm(i) = 0 - mp2Elm(i) = INVALID_ELM_ID - EXIT - endif - end do - - if (inBound) then - xc = 0.0_MPAS_RKIND - yc = 0.0_MPAS_RKIND - zc = 0.0_MPAS_RKIND - do k = 1, nEdgesOnCell(i) - j = verticesOnCell(k,i) - xc = xc + xVertex(j) - yc = yc + yVertex(j) - zc = zc + zVertex(j) - xComputed = sphereRadius*cos(latVertex(j))*cos(lonVertex(j)) - yComputed = sphereRadius*cos(latVertex(j))*sin(lonVertex(j)) - zComputed = sphereRadius*sin(latVertex(j)) - latComputed = asin(zVertex(j)/sphereRadius) - lonComputed = atan2(yVertex(j),xVertex(j)) - if (lonComputed .le. 0.0_MPAS_RKIND) then ! lon[0,2pi] - lonComputed = lonComputed + 2.0_MPAS_RKIND*pi - endif - - end do - xc = xc/nEdgesOnCell(i) - yc = yc/nEdgesOnCell(i) - zc = zc/nEdgesOnCell(i) - ! normalize - radius = sqrt(xc*xc + yc*yc + zc*zc)! assuming sphere center to be at origin - xc = xc/radius * sphereRadius - yc = yc/radius * sphereRadius - zc = zc/radius * sphereRadius - mpPosition(1,i) = xc - mpPosition(2,i) = yc - mpPosition(3,i) = zc - mpLatLon(1,i) = asin(zc/sphereRadius) - lon = atan2(yc,xc) - if (lon .le. 0.0_MPAS_RKIND) then ! lon[0,2pi] - lon = lon + 2.0_MPAS_RKIND*pi - endif - mpLatLon(2,i) = lon - endif - end do - ! check first element/cell for delta - maxlon = minval(lonVertex) - minlon = maxval(lonVertex) - do i = 1, nEdgesOnCell(1) - j = verticesOnCell(i,1) - if(maxlon .lt. lonVertex(j)) then - maxlon = lonVertex(j) - endif - if(minlon .gt. lonVertex(j)) then - minlon = lonVertex(j) - endif - end do - 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)) - - - deltaLon = maxlon - minlon - - do i = 1,nVertices - dispIncr(1,i) = sphereRadius*cos(latVertex(i))*deltaLon - dispIncr(2,i) = 0.0_MPAS_RKIND - end do - call polympo_setMeshOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr)) - call polympo_push(mpMesh) - do i = 1,nVertices - dispIncr(1,i) = sphereRadius*cos(latVertex(i))*2*deltaLon - dispIncr(2,i) = 0.0_MPAS_RKIND - end do - call polympo_setMeshOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr)) - call polympo_push(mpMesh) - call polympo_deleteMPMesh(mpMesh) - call polympo_finalize() - - call mpi_finalize(ierr) - - deallocate(nEdgesOnCell) - deallocate(xVertex) - deallocate(yVertex) - deallocate(zVertex) - deallocate(latVertex) - deallocate(lonVertex) - deallocate(verticesOnCell) - deallocate(cellsOnCell) - deallocate(dispIncr) - deallocate(mpsPerElm) - deallocate(mp2Elm) - deallocate(isMPActive) - deallocate(mpPosition) - deallocate(mpLatLon) + ! call mpi_init(ierr) + ! call mpi_comm_rank(mpi_comm_handle, self, ierr) + + ! call polympo_setMPICommunicator(mpi_comm_handle) + ! call polympo_initialize() + ! call polympo_enableTiming() + + ! call polympo_checkPrecisionForRealKind(MPAS_RKIND) + ! argc = command_argument_count() + ! if(argc == 3) then + ! call get_command_argument(1, input) + ! read(input, '(I7)') mpsScaleFactorPerVtx + ! call get_command_argument(2, input) + ! read(input, '(I7)') numPush + ! call get_command_argument(3, filename) + ! else + ! write(0, *) "Usage: ./testFortranMPAdvection " + ! end if + + ! call readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & + ! nCells, nVertices, nEdgesOnCell, & + ! onSphere, sphereRadius, & + ! xVertex, yVertex, zVertex, & + ! latVertex, lonVertex, & + ! verticesOnCell, cellsOnCell) + ! if (onSphere .ne. 'YES') then + ! write (*,*) "The mesh is not spherical!" + ! call exit(1) + ! end if + + ! setMeshOption = 0 !create an empty mesh + ! setMPOption = 0 !create an empty set of MPs + ! mpMesh = polympo_createMPMesh(setMeshOption,setMPOption) !creates test mesh + ! call loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & + ! nCells, nVertices, nEdgesOnCell, & + ! onSphere, sphereRadius, & + ! xVertex, yVertex, zVertex, & + ! latVertex, & + ! verticesOnCell, cellsOnCell) + + ! !createMPs + ! numMPs = 0 + ! do i = 1, nCells + ! numMPs = numMPs + nEdgesOnCell(i) * mpsScaleFactorPerVtx + ! end do + + ! print *, "Scale Factor", mpsScaleFactorPerVtx + ! print *, "NUM MPs", numMPs + + ! allocate(mpsPerElm(nCells)) + ! allocate(mp2Elm(numMPs)) + ! allocate(isMPActive(numMPs)) + ! allocate(mpPosition(3,numMPs)) + ! allocate(mpLatLon(2,numMPs)) + + ! isMPActive = MP_ACTIVE !all active MPs and some changed below + + ! numMPsCount = 0 + ! do i = 1, nCells + ! localNumMPs = nEdgesOnCell(i) * mpsScaleFactorPerVtx + ! mp2Elm(numMPsCount+1:numMPsCount+localNumMPs) = i + ! mpsPerElm(i) = localNumMPs + ! numMPsCount = numMPsCount + localNumMPs + ! end do + + ! call assert(numMPsCount == numMPs, "num mps miscounted") + + ! numMPsCount = 0 + ! do i = 1, nCells + ! xc = 0.0_MPAS_RKIND + ! yc = 0.0_MPAS_RKIND + ! zc = 0.0_MPAS_RKIND + ! do k = 1, nEdgesOnCell(i) + ! j = verticesOnCell(k,i) + ! xc = xc + xVertex(j) + ! yc = yc + yVertex(j) + ! zc = zc + zVertex(j) + ! end do + ! xc = xc/nEdgesOnCell(i) + ! yc = yc/nEdgesOnCell(i) + ! zc = zc/nEdgesOnCell(i) + + ! do k = 1, nEdgesOnCell(i) + ! j = verticesOnCell(k,i) + + ! ! note: m=0 not included but should lead to x_mp=xc and m=mpsScaleFactorPerVtx+1 is also not include but should lead to x_mp=xVertex(j) + ! ! so'mpsScaleFactorPerVtx+1' segments and 'mpsScaleFactorPerVtx+2' points along line from 'xc' to 'xVertex(j)' + ! ! taking only "inner" or interior 'mpsScaleFactorPerVtx' points (i.e., exclude end points of 'xc' and 'xVertex(j)') and same applies to y- and z-coordinates + ! do m = 1, mpsScaleFactorPerVtx + ! numMPsCount = numMPsCount + 1 + ! xMP = (mpsScaleFactorPerVtx+1 - m) * xc + m*xVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation + ! yMP = (mpsScaleFactorPerVtx+1 - m) * yc + m*yVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation + ! zMP = (mpsScaleFactorPerVtx+1 - m) * zc + m*zVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation + + ! ! normalize to project each MP to be on sphere of radius 'sphereRadius' + ! radius = sqrt(xMP*xMP + yMP*yMP + zMP*zMP) ! assuming sphere center to be at origin + ! xMP = xMP/radius * sphereRadius + ! yMP = yMP/radius * sphereRadius + ! zMP = zMP/radius * sphereRadius + ! mpPosition(1,numMPsCount) = xMP + ! mpPosition(2,numMPsCount) = yMP + ! mpPosition(3,numMPsCount) = zMP + ! mpLatLon(1,numMPsCount) = asin(zMP/sphereRadius) + ! lon = atan2(yMP,xMP) + ! if (lon .le. 0.0_MPAS_RKIND) then ! lon[0,2pi] + ! lon = lon + 2.0_MPAS_RKIND*pi + ! endif + ! mpLatLon(2,numMPsCount) = lon + ! end do + ! end do + ! end do + + ! call assert(numMPsCount == numMPs, "num mps miscounted") + + ! 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)) + + ! call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) + + ! call polympo_summarizeTime(); + + ! call polympo_deleteMPMesh(mpMesh) + ! call polympo_finalize() + + ! call mpi_finalize(ierr) + + ! deallocate(nEdgesOnCell) + ! deallocate(xVertex) + ! deallocate(yVertex) + ! deallocate(zVertex) + ! deallocate(latVertex) + ! deallocate(lonVertex) + ! deallocate(verticesOnCell) + ! deallocate(cellsOnCell) + ! deallocate(mpsPerElm) + ! deallocate(mp2Elm) + ! deallocate(isMPActive) + ! deallocate(mpPosition) + ! deallocate(mpLatLon) stop + end program From 0232f3c7955a14e32154aeacb94354cdac5c39b1 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Mon, 15 Jul 2024 20:21:14 -0700 Subject: [PATCH 15/26] calcsurfDispIncr --- src/pmpo_c.cpp | 2 +- src/pmpo_c.h | 2 +- src/pmpo_fortran.f90 | 4 +- test/calculateDisplacement.f90 | 68 ++++++++++++++++---------------- test/testFortran.f90 | 2 +- test/testFortranMPMeshModule.f90 | 2 +- 6 files changed, 40 insertions(+), 40 deletions(-) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index c447dd6d..a71922c6 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -595,7 +595,7 @@ void polympo_getMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, cons Kokkos::deep_copy(arrayHost, vtxField); } -void polympo_setMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) { +void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) { //check mpMesh is valid checkMPMeshValid(p_mpmesh); auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; diff --git a/src/pmpo_c.h b/src/pmpo_c.h index dfd6116b..5390dcda 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -59,7 +59,7 @@ void polympo_setMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, const void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double* latitude); void polympo_setMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d void polympo_getMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d -void polympo_setMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d +void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d void polympo_getMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d // calculations diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index d7d7b9a3..761bbb44 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -483,8 +483,8 @@ subroutine polympo_getMeshOnSurfVeloIncr(mpMesh, nComps, nVertices, array) & !> @param nVertices(in) numVertices !> @param array(in) input mesh velocity 2D array (2,numVtx) !--------------------------------------------------------------------------- - subroutine polympo_setMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) & - bind(C, NAME='polympo_setMeshOnSurfDispIncr_f') + subroutine polympo_setMeshVtxOnSurfDispIncr(mpMesh, nComps, nVertices, array) & + bind(C, NAME='polympo_setMeshVtxOnSurfDispIncr_f') use :: iso_c_binding type(c_ptr), value :: mpMesh integer(c_int), value :: nComps, nVertices diff --git a/test/calculateDisplacement.f90 b/test/calculateDisplacement.f90 index dac64c2c..c499008b 100644 --- a/test/calculateDisplacement.f90 +++ b/test/calculateDisplacement.f90 @@ -1,40 +1,40 @@ -! subroutine calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) -! use :: polympo -! use :: readMPAS -! use :: iso_c_binding -! implicit none +subroutine calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) + use :: polympo + use :: readMPAS + use :: iso_c_binding + implicit none -! real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex -! real(kind=MPAS_RKIND) :: maxlon, minlon, deltaLon, sphereRadius -! integer, dimension(:), pointer :: nEdgesOnCell -! integer, dimension(:,:), pointer :: verticesOnCell -! integer :: i, j, nVertices, nCompsDisp -! real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr -! type(c_ptr) :: mpMesh + real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND) :: maxlon, minlon, deltaLon, sphereRadius + integer, dimension(:), pointer :: nEdgesOnCell + integer, dimension(:,:), pointer :: verticesOnCell + integer :: i, j, nVertices, nCompsDisp + real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr + type(c_ptr) :: mpMesh -! nCompsDisp = 2 -! allocate(dispIncr(nCompsDisp,nVertices)) + nCompsDisp = 2 + allocate(dispIncr(nCompsDisp,nVertices)) -! ! check first element/cell for delta -! maxlon = minval(lonVertex) -! minlon = maxval(lonVertex) -! do i = 1, nEdgesOnCell(1) -! j = verticesOnCell(i,1) -! if(maxlon .lt. lonVertex(j)) then -! maxlon = lonVertex(j) -! endif -! if(minlon .gt. lonVertex(j)) then -! minlon = lonVertex(j) -! endif -! end do + ! check first element/cell for delta + maxlon = minval(lonVertex) + minlon = maxval(lonVertex) + do i = 1, nEdgesOnCell(1) + j = verticesOnCell(i,1) + if(maxlon .lt. lonVertex(j)) then + maxlon = lonVertex(j) + endif + if(minlon .gt. lonVertex(j)) then + minlon = lonVertex(j) + endif + end do -! deltaLon = maxlon - minlon + deltaLon = maxlon - minlon -! do i = 1,nVertices -! dispIncr(1,i) = sphereRadius*cos(latVertex(i))*deltaLon -! dispIncr(2,i) = 0.0_MPAS_RKIND -! end do -! call polympo_setMeshVtxOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr)) + do i = 1,nVertices + dispIncr(1,i) = sphereRadius*cos(latVertex(i))*deltaLon + dispIncr(2,i) = 0.0_MPAS_RKIND + end do + call polympo_setMeshVtxOnSurfDispIncr(mpMesh,nCompsDisp,nVertices,c_loc(dispIncr)) -! deallocate(dispIncr) -! end subroutine \ No newline at end of file + deallocate(dispIncr) +end subroutine \ No newline at end of file diff --git a/test/testFortran.f90 b/test/testFortran.f90 index dd953ace..92ecb2ea 100644 --- a/test/testFortran.f90 +++ b/test/testFortran.f90 @@ -75,7 +75,7 @@ program main end do end do call polympo_setMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray)) - call polympo_setMeshOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray)) + call polympo_setMeshVtxOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray)) Mesharray = 1 call polympo_getMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray)) diff --git a/test/testFortranMPMeshModule.f90 b/test/testFortranMPMeshModule.f90 index 6c90fdf0..ee02098b 100644 --- a/test/testFortranMPMeshModule.f90 +++ b/test/testFortranMPMeshModule.f90 @@ -76,7 +76,7 @@ program main Mesharray = value1 call polympo_setMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray)) Mesharray = value2 - call polympo_setMeshOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray)) + call polympo_setMeshVtxOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray)) Mesharray = 1 call polympo_getMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray)) From 6f11255fa20efce8af4258b42ad9a0bd4e9335b9 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Mon, 15 Jul 2024 20:35:44 -0700 Subject: [PATCH 16/26] advection passing --- test/testFortranMPAdvection.f90 | 272 ++++++++++++++++---------------- 1 file changed, 136 insertions(+), 136 deletions(-) diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 6f7755fa..eefc4747 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -16,10 +16,10 @@ subroutine runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, real(kind=MPAS_RKIND) :: sphereRadius - ! do i = 1, numPush - ! call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) - ! call polympo_push(mpMesh) - ! end do + do i = 1, numPush + call calcSurfDispIncr(mpMesh, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) + call polympo_push(mpMesh) + end do end subroutine @@ -58,146 +58,146 @@ program main integer, parameter :: MP_INACTIVE = 0 integer, parameter :: INVALID_ELM_ID = -1 - ! call mpi_init(ierr) - ! call mpi_comm_rank(mpi_comm_handle, self, ierr) + call mpi_init(ierr) + call mpi_comm_rank(mpi_comm_handle, self, ierr) - ! call polympo_setMPICommunicator(mpi_comm_handle) - ! call polympo_initialize() + call polympo_setMPICommunicator(mpi_comm_handle) + call polympo_initialize() ! call polympo_enableTiming() - ! call polympo_checkPrecisionForRealKind(MPAS_RKIND) - ! argc = command_argument_count() - ! if(argc == 3) then - ! call get_command_argument(1, input) - ! read(input, '(I7)') mpsScaleFactorPerVtx - ! call get_command_argument(2, input) - ! read(input, '(I7)') numPush - ! call get_command_argument(3, filename) - ! else - ! write(0, *) "Usage: ./testFortranMPAdvection " - ! end if - - ! call readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & - ! nCells, nVertices, nEdgesOnCell, & - ! onSphere, sphereRadius, & - ! xVertex, yVertex, zVertex, & - ! latVertex, lonVertex, & - ! verticesOnCell, cellsOnCell) - ! if (onSphere .ne. 'YES') then - ! write (*,*) "The mesh is not spherical!" - ! call exit(1) - ! end if - - ! setMeshOption = 0 !create an empty mesh - ! setMPOption = 0 !create an empty set of MPs - ! mpMesh = polympo_createMPMesh(setMeshOption,setMPOption) !creates test mesh - ! call loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & - ! nCells, nVertices, nEdgesOnCell, & - ! onSphere, sphereRadius, & - ! xVertex, yVertex, zVertex, & - ! latVertex, & - ! verticesOnCell, cellsOnCell) - - ! !createMPs - ! numMPs = 0 - ! do i = 1, nCells - ! numMPs = numMPs + nEdgesOnCell(i) * mpsScaleFactorPerVtx - ! end do - - ! print *, "Scale Factor", mpsScaleFactorPerVtx - ! print *, "NUM MPs", numMPs - - ! allocate(mpsPerElm(nCells)) - ! allocate(mp2Elm(numMPs)) - ! allocate(isMPActive(numMPs)) - ! allocate(mpPosition(3,numMPs)) - ! allocate(mpLatLon(2,numMPs)) - - ! isMPActive = MP_ACTIVE !all active MPs and some changed below - - ! numMPsCount = 0 - ! do i = 1, nCells - ! localNumMPs = nEdgesOnCell(i) * mpsScaleFactorPerVtx - ! mp2Elm(numMPsCount+1:numMPsCount+localNumMPs) = i - ! mpsPerElm(i) = localNumMPs - ! numMPsCount = numMPsCount + localNumMPs - ! end do - - ! call assert(numMPsCount == numMPs, "num mps miscounted") - - ! numMPsCount = 0 - ! do i = 1, nCells - ! xc = 0.0_MPAS_RKIND - ! yc = 0.0_MPAS_RKIND - ! zc = 0.0_MPAS_RKIND - ! do k = 1, nEdgesOnCell(i) - ! j = verticesOnCell(k,i) - ! xc = xc + xVertex(j) - ! yc = yc + yVertex(j) - ! zc = zc + zVertex(j) - ! end do - ! xc = xc/nEdgesOnCell(i) - ! yc = yc/nEdgesOnCell(i) - ! zc = zc/nEdgesOnCell(i) - - ! do k = 1, nEdgesOnCell(i) - ! j = verticesOnCell(k,i) + call polympo_checkPrecisionForRealKind(MPAS_RKIND) + argc = command_argument_count() + if(argc == 3) then + call get_command_argument(1, input) + read(input, '(I7)') mpsScaleFactorPerVtx + call get_command_argument(2, input) + read(input, '(I7)') numPush + call get_command_argument(3, filename) + else + write(0, *) "Usage: ./testFortranMPAdvection " + end if + + call readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & + nCells, nVertices, nEdgesOnCell, & + onSphere, sphereRadius, & + xVertex, yVertex, zVertex, & + latVertex, lonVertex, & + verticesOnCell, cellsOnCell) + if (onSphere .ne. 'YES') then + write (*,*) "The mesh is not spherical!" + call exit(1) + end if + + setMeshOption = 0 !create an empty mesh + setMPOption = 0 !create an empty set of MPs + mpMesh = polympo_createMPMesh(setMeshOption,setMPOption) !creates test mesh + call loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, & + nCells, nVertices, nEdgesOnCell, & + onSphere, sphereRadius, & + xVertex, yVertex, zVertex, & + latVertex, & + verticesOnCell, cellsOnCell) + + !createMPs + numMPs = 0 + do i = 1, nCells + numMPs = numMPs + nEdgesOnCell(i) * mpsScaleFactorPerVtx + end do + + print *, "Scale Factor", mpsScaleFactorPerVtx + print *, "NUM MPs", numMPs + + allocate(mpsPerElm(nCells)) + allocate(mp2Elm(numMPs)) + allocate(isMPActive(numMPs)) + allocate(mpPosition(3,numMPs)) + allocate(mpLatLon(2,numMPs)) + + isMPActive = MP_ACTIVE !all active MPs and some changed below + + numMPsCount = 0 + do i = 1, nCells + localNumMPs = nEdgesOnCell(i) * mpsScaleFactorPerVtx + mp2Elm(numMPsCount+1:numMPsCount+localNumMPs) = i + mpsPerElm(i) = localNumMPs + numMPsCount = numMPsCount + localNumMPs + end do + + call assert(numMPsCount == numMPs, "num mps miscounted") + + numMPsCount = 0 + do i = 1, nCells + xc = 0.0_MPAS_RKIND + yc = 0.0_MPAS_RKIND + zc = 0.0_MPAS_RKIND + do k = 1, nEdgesOnCell(i) + j = verticesOnCell(k,i) + xc = xc + xVertex(j) + yc = yc + yVertex(j) + zc = zc + zVertex(j) + end do + xc = xc/nEdgesOnCell(i) + yc = yc/nEdgesOnCell(i) + zc = zc/nEdgesOnCell(i) + + do k = 1, nEdgesOnCell(i) + j = verticesOnCell(k,i) - ! ! note: m=0 not included but should lead to x_mp=xc and m=mpsScaleFactorPerVtx+1 is also not include but should lead to x_mp=xVertex(j) - ! ! so'mpsScaleFactorPerVtx+1' segments and 'mpsScaleFactorPerVtx+2' points along line from 'xc' to 'xVertex(j)' - ! ! taking only "inner" or interior 'mpsScaleFactorPerVtx' points (i.e., exclude end points of 'xc' and 'xVertex(j)') and same applies to y- and z-coordinates - ! do m = 1, mpsScaleFactorPerVtx - ! numMPsCount = numMPsCount + 1 - ! xMP = (mpsScaleFactorPerVtx+1 - m) * xc + m*xVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation - ! yMP = (mpsScaleFactorPerVtx+1 - m) * yc + m*yVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation - ! zMP = (mpsScaleFactorPerVtx+1 - m) * zc + m*zVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation + ! note: m=0 not included but should lead to x_mp=xc and m=mpsScaleFactorPerVtx+1 is also not include but should lead to x_mp=xVertex(j) + ! so'mpsScaleFactorPerVtx+1' segments and 'mpsScaleFactorPerVtx+2' points along line from 'xc' to 'xVertex(j)' + ! taking only "inner" or interior 'mpsScaleFactorPerVtx' points (i.e., exclude end points of 'xc' and 'xVertex(j)') and same applies to y- and z-coordinates + do m = 1, mpsScaleFactorPerVtx + numMPsCount = numMPsCount + 1 + xMP = (mpsScaleFactorPerVtx+1 - m) * xc + m*xVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation + yMP = (mpsScaleFactorPerVtx+1 - m) * yc + m*yVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation + zMP = (mpsScaleFactorPerVtx+1 - m) * zc + m*zVertex(j) / (mpsScaleFactorPerVtx+1) ! linear interpolation - ! ! normalize to project each MP to be on sphere of radius 'sphereRadius' - ! radius = sqrt(xMP*xMP + yMP*yMP + zMP*zMP) ! assuming sphere center to be at origin - ! xMP = xMP/radius * sphereRadius - ! yMP = yMP/radius * sphereRadius - ! zMP = zMP/radius * sphereRadius - ! mpPosition(1,numMPsCount) = xMP - ! mpPosition(2,numMPsCount) = yMP - ! mpPosition(3,numMPsCount) = zMP - ! mpLatLon(1,numMPsCount) = asin(zMP/sphereRadius) - ! lon = atan2(yMP,xMP) - ! if (lon .le. 0.0_MPAS_RKIND) then ! lon[0,2pi] - ! lon = lon + 2.0_MPAS_RKIND*pi - ! endif - ! mpLatLon(2,numMPsCount) = lon - ! end do - ! end do - ! end do - - ! call assert(numMPsCount == numMPs, "num mps miscounted") - - ! 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)) - - ! call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) + ! normalize to project each MP to be on sphere of radius 'sphereRadius' + radius = sqrt(xMP*xMP + yMP*yMP + zMP*zMP) ! assuming sphere center to be at origin + xMP = xMP/radius * sphereRadius + yMP = yMP/radius * sphereRadius + zMP = zMP/radius * sphereRadius + mpPosition(1,numMPsCount) = xMP + mpPosition(2,numMPsCount) = yMP + mpPosition(3,numMPsCount) = zMP + mpLatLon(1,numMPsCount) = asin(zMP/sphereRadius) + lon = atan2(yMP,xMP) + if (lon .le. 0.0_MPAS_RKIND) then ! lon[0,2pi] + lon = lon + 2.0_MPAS_RKIND*pi + endif + mpLatLon(2,numMPsCount) = lon + end do + end do + end do + + call assert(numMPsCount == numMPs, "num mps miscounted") + + 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)) + + call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) ! call polympo_summarizeTime(); - ! call polympo_deleteMPMesh(mpMesh) - ! call polympo_finalize() - - ! call mpi_finalize(ierr) - - ! deallocate(nEdgesOnCell) - ! deallocate(xVertex) - ! deallocate(yVertex) - ! deallocate(zVertex) - ! deallocate(latVertex) - ! deallocate(lonVertex) - ! deallocate(verticesOnCell) - ! deallocate(cellsOnCell) - ! deallocate(mpsPerElm) - ! deallocate(mp2Elm) - ! deallocate(isMPActive) - ! deallocate(mpPosition) - ! deallocate(mpLatLon) + call polympo_deleteMPMesh(mpMesh) + call polympo_finalize() + + call mpi_finalize(ierr) + + deallocate(nEdgesOnCell) + deallocate(xVertex) + deallocate(yVertex) + deallocate(zVertex) + deallocate(latVertex) + deallocate(lonVertex) + deallocate(verticesOnCell) + deallocate(cellsOnCell) + deallocate(mpsPerElm) + deallocate(mp2Elm) + deallocate(isMPActive) + deallocate(mpPosition) + deallocate(mpLatLon) stop From 5d062e18f95a3c3fd5899d0d9e486b8cb2d2e389 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Tue, 16 Jul 2024 12:27:19 -0700 Subject: [PATCH 17/26] pumipic timing --- src/pmpo_c.cpp | 12 ++++++++++++ src/pmpo_c.h | 5 +++++ src/pmpo_fortran.f90 | 19 +++++++++++++++++++ test/CMakeLists.txt | 2 +- test/testFortranMPAdvection.f90 | 5 +++-- 5 files changed, 40 insertions(+), 3 deletions(-) diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp index a71922c6..4991c80e 100644 --- a/src/pmpo_c.cpp +++ b/src/pmpo_c.cpp @@ -655,4 +655,16 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a Kokkos::View owningProc("owningProc",nCells); Kokkos::deep_copy(owningProc, arrayHost); p_mesh->setOwningProc(owningProc); +} + +void polympo_enableTiming_f(){ + pumipic::EnableTiming(); +} + +void polympo_summarizeTime_f(){ + pumipic::SummarizeTime(); +} + +void polympo_setTimingVerbosity_f(int v){ + pumipic::SetTimingVerbosity(v); } \ No newline at end of file diff --git a/src/pmpo_c.h b/src/pmpo_c.h index 5390dcda..db8e35cd 100644 --- a/src/pmpo_c.h +++ b/src/pmpo_c.h @@ -64,5 +64,10 @@ void polympo_getMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, cons // calculations void polympo_push_f(MPMesh_ptr p_mpmesh); + +// Timing +void polympo_enableTiming_f(); +void polympo_summarizeTime_f(); +void polympo_setTimingVerbosity_f(int v); } #endif diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90 index 761bbb44..4753e179 100644 --- a/src/pmpo_fortran.f90 +++ b/src/pmpo_fortran.f90 @@ -530,6 +530,25 @@ subroutine polympo_push(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 + !--------------------------------------------------------------------------- + subroutine polympo_setTimingVerbosity(v) bind(C, NAME='polympo_setTimingVerbosity_f') + use :: iso_c_binding + integer(c_int), value :: v + end subroutine end interface contains !--------------------------------------------------------------------------- diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 533a0055..7376e265 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -87,7 +87,7 @@ if(EXISTS ${TEST_NC_FILE_SPHERICAL}) message(STATUS "Found Spherical NC File: ${TEST_NC_FILE_SPHERICAL}") set(TEST_NC_FILE ${TEST_NC_FILE_SPHERICAL}) # this test only works for spherical NC files! - pmpo_add_test(testFortranMPAdvection 1 ./testFortranMPAdvection ${TEST_NC_FILE}) + pmpo_add_test(testFortranMPAdvection 1 ./testFortranMPAdvection 5 5 ${TEST_NC_FILE}) else() set(TEST_NC_FILE ${TEST_NC_FILE_PLANAR}) endif() diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index eefc4747..7798f3c0 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -63,7 +63,7 @@ program main call polympo_setMPICommunicator(mpi_comm_handle) call polympo_initialize() - ! call polympo_enableTiming() + call polympo_enableTiming() call polympo_checkPrecisionForRealKind(MPAS_RKIND) argc = command_argument_count() @@ -75,6 +75,7 @@ program main call get_command_argument(3, filename) else write(0, *) "Usage: ./testFortranMPAdvection " + call exit(1) end if call readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, & @@ -178,7 +179,7 @@ program main call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) - ! call polympo_summarizeTime(); + call polympo_summarizeTime(); call polympo_deleteMPMesh(mpMesh) call polympo_finalize() From 6ab0ae80c677405dc2bd4a9bfe35581bd5d88e7c Mon Sep 17 00:00:00 2001 From: Angelyr Date: Tue, 16 Jul 2024 16:11:02 -0700 Subject: [PATCH 18/26] set Proc Wedges --- test/testFortranMPAdvection.f90 | 39 +++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 7798f3c0..9cbf2ec1 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -2,6 +2,45 @@ module advectionTests contains include "calculateDisplacement.f90" + subroutine setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, verticesOnCell) + use :: polympo + use :: readMPAS + use :: iso_c_binding + implicit none + + type(c_ptr) :: mpMesh + integer :: i, j, k, nVertices, nCells, comm_size + integer, dimension(:), pointer :: owningProc, nEdgesOnCell + integer, dimension(:,:), pointer :: verticesOnCell + real(kind=MPAS_RKIND), dimension(:), pointer :: elmLat, vtxRotLat + real(kind=MPAS_RKIND) :: normalizedLat, sum, min, max + + allocate(elmLat(nCells)) + allocate(owningProc(nCells)) + allocate(vtxRotLat(nVertices)) + + call polympo_getMeshVtxRotLat(mpMesh, nVertices, c_loc(vtxRotLat)) + + min = 0 + max = 0 + do i = 1, nCells + sum = 0; + do j = 1, nEdgesOnCell(i) + sum = sum + vtxRotLat(verticesOnCell(j,i)) + end do + elmLat(i) = sum/nEdgesOnCell(i); + if (elmLat(i) < min) min = elmLat(i) + if (elmLat(i) > max) max = elmLat(i) + end do + + do i = 1, nCells + normalizedLat = (elmLat(i) - min) / (max - min); + owningProc(i) = normalizedLat * comm_size; + end do + + ! TODO SET OWNING PROC + end subroutine + subroutine runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) use :: polympo use :: readMPAS From bb00e558170ade14f9c59cb4b72f5bd71c38a01e Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 17 Jul 2024 11:50:22 -0700 Subject: [PATCH 19/26] fixed mpi test perlmutter --- test/CMakeLists.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7376e265..500fb58e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -58,12 +58,12 @@ pmpo_add_fortran_exe(testFortranMPAdvection testFortranMPAdvection.f90) pmpo_add_fortran_exe(testFortranCreateRebuildMPs testFortranCreateRebuildMPs.f90) #add test -find_program(MPIRUN_EXECUTABLE NAMES mpirun) +find_program(MPIRUN NAMES mpirun mpiexec srun) function(pmpo_add_test testname procs command) - if(MPIRUN_EXECUTABLE) + if(MPIRUN) add_test(NAME ${testname} - COMMAND mpirun --bind-to core -np ${procs} ${command} ${ARGN}) + COMMAND ${MPIRUN} -n ${procs} ${command} ${ARGN}) else() add_test(NAME ${testname} COMMAND ${command} ${ARGN}) endif() @@ -96,6 +96,6 @@ pmpo_add_test(testFortranCreateRebuildMPs 1 ./testFortranCreateRebuildMPs ${TEST pmpo_add_test(test_timing 1 ./timeAssmblyWachspress 1 ${TEST_NC_FILE}) pmpo_add_test(test_wachspress 1 ./testWachspress ${TEST_NC_FILE}) -if(MPIRUN_EXECUTABLE) +if(MPIRUN) pmpo_add_test(test_migration 8 ./testMigration ${TEST_NC_FILE}) endif() \ No newline at end of file From 6bf800318c735ad87020f26615cb5451575a71c4 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 17 Jul 2024 12:42:36 -0700 Subject: [PATCH 20/26] set owning proc --- test/testFortranMPAdvection.f90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 9cbf2ec1..ed132195 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -38,7 +38,8 @@ subroutine setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, ver owningProc(i) = normalizedLat * comm_size; end do - ! TODO SET OWNING PROC + call polympo_startMeshFill(mpMesh) + call polympo_setOwningProc(mpMesh, nCells, c_loc(owningProc)) end subroutine subroutine runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) @@ -76,7 +77,7 @@ program main !integer, parameter :: APP_RKIND = selected_real_kind(15) type(c_ptr) :: mpMesh - integer :: ierr, self + integer :: ierr, self, comm_size integer :: argc, i, j, arglen, k, m, mpsScaleFactorPerVtx, localNumMPs integer :: setMeshOption, setMPOption integer :: maxEdges, vertexDegree, nCells, nVertices @@ -99,6 +100,7 @@ program main call mpi_init(ierr) call mpi_comm_rank(mpi_comm_handle, self, ierr) + call mpi_comm_size(mpi_comm_handle, comm_size, ierr) call polympo_setMPICommunicator(mpi_comm_handle) call polympo_initialize() @@ -216,6 +218,7 @@ program main call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon)) call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition)) + call setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, verticesOnCell) call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) call polympo_summarizeTime(); From 4448cef7276a4d00f480bd018caa2d3cb4646c02 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 17 Jul 2024 13:01:23 -0700 Subject: [PATCH 21/26] fixed owning proc --- test/testFortranMPAdvection.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index ed132195..476f9830 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -34,8 +34,9 @@ subroutine setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, ver end do do i = 1, nCells - normalizedLat = (elmLat(i) - min) / (max - min); - owningProc(i) = normalizedLat * comm_size; + normalizedLat = (elmLat(i) - min) / (max - min) + owningProc(i) = normalizedLat * comm_size + owningProc(i) = 0 end do call polympo_startMeshFill(mpMesh) From 9733e997b7390569f2bb50db744997b0f7c2d756 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 17 Jul 2024 13:33:30 -0700 Subject: [PATCH 22/26] timers --- src/pmpo_MPMesh.cpp | 6 ++++++ src/pmpo_materialPoints.cpp | 2 +- src/pmpo_materialPoints.hpp | 2 ++ src/pmpo_mesh.cpp | 2 ++ src/pmpo_mesh.hpp | 1 + src/pmpo_wachspressBasis.hpp | 2 ++ 6 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index eed994a1..5922f11d 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -74,6 +74,7 @@ void MPMesh::CVTTrackingEdgeCenterBased(Vec2dView dx){ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ + Kokkos::Timer timer; int numVtxs = p_mesh->getNumVertices(); int numElms = p_mesh->getNumElements(); auto numMPs = p_MPs->getCount(); @@ -197,6 +198,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){ fprintf(pFile," \n \n \n \n\n"); fclose(pFile); } + pumipic::RecordTime("PolyMPO_CVTTrackingElmCenterBased", timer.seconds()); } void MPMesh::T2LTracking(Vec2dView dx){ @@ -263,6 +265,7 @@ void MPMesh::T2LTracking(Vec2dView dx){ } bool getAnyIsMigrating(bool isMigrating) { + Kokkos::Timer timer; int comm_rank; MPI_Comm_rank(MPI_COMM_WORLD, &comm_rank); int comm_size; @@ -270,10 +273,12 @@ bool getAnyIsMigrating(bool isMigrating) { bool anyIsMigrating = false; MPI_Allreduce(&isMigrating, &anyIsMigrating, 1, MPI_C_BOOL, MPI_LOR, MPI_COMM_WORLD); + pumipic::RecordTime("PolyMPO_getAnyIsMigrating", timer.seconds()); return anyIsMigrating; } void MPMesh::push(){ + Kokkos::Timer timer; p_mesh->computeRotLatLonIncr(); sphericalInterpolation(*this); p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ @@ -289,6 +294,7 @@ void MPMesh::push(){ p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ p_MPs->updateMPSlice(); // Tgt becomes Cur + pumipic::RecordTime("PolyMPO_push", timer.seconds()); } void MPMesh::printVTP_mesh(int printVTPIndex){ diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp index a9430dc8..c1570ef8 100644 --- a/src/pmpo_materialPoints.cpp +++ b/src/pmpo_materialPoints.cpp @@ -124,7 +124,7 @@ bool MaterialPoints::migrate() { if (getOpMode() == polyMPO::MP_DEBUG) printf("Material point migration: %f\n", timer.seconds()); - + pumipic::RecordTime("PolyMPO_migrate", timer.seconds()); return pumipic::getLastValue(isMigrating) > 0; } diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp index 29e56ee1..ae93a2e4 100644 --- a/src/pmpo_materialPoints.hpp +++ b/src/pmpo_materialPoints.hpp @@ -208,6 +208,7 @@ class MaterialPoints { updateMPSlice(); } void updateRotLatLonAndXYZ2Tgt(const double radius){ + Kokkos::Timer timer; auto curPosRotLatLon = MPs->get(); auto tgtPosRotLatLon = MPs->get(); auto tgtPosXYZ = MPs->get(); @@ -233,6 +234,7 @@ class MaterialPoints { PMT_ALWAYS_ASSERT(false); } ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude"); + pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds()); } template diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp index 084e4352..269aa4cb 100644 --- a/src/pmpo_mesh.cpp +++ b/src/pmpo_mesh.cpp @@ -40,6 +40,7 @@ namespace polyMPO{ } void Mesh::computeRotLatLonIncr(){ + Kokkos::Timer timer; PMT_ALWAYS_ASSERT(geomType_ == geom_spherical_surf); auto dispIncr = getMeshField(); @@ -53,6 +54,7 @@ namespace polyMPO{ 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() { diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp index 22e4d8b2..9e82c741 100644 --- a/src/pmpo_mesh.hpp +++ b/src/pmpo_mesh.hpp @@ -2,6 +2,7 @@ #define POLYMPO_MESH_H #include "pmpo_utils.hpp" +#include namespace polyMPO{ diff --git a/src/pmpo_wachspressBasis.hpp b/src/pmpo_wachspressBasis.hpp index cff2a35a..a4c64652 100644 --- a/src/pmpo_wachspressBasis.hpp +++ b/src/pmpo_wachspressBasis.hpp @@ -306,6 +306,7 @@ 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(); @@ -355,6 +356,7 @@ void sphericalInterpolation(MPMesh& mpMesh){ } }; p_MPs->parallel_for(interpolation, "interpolation"); + pumipic::RecordTime("PolyMPO_Reconstruct", timer.seconds()); } } //namespace polyMPO end From 96760034b69be429255cd89551661b5758973cac Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 17 Jul 2024 19:46:38 -0700 Subject: [PATCH 23/26] fixed error --- test/testFortranMPAdvection.f90 | 6 +++--- test/testMigration.cpp | 1 + 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 476f9830..61ee98db 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -21,8 +21,8 @@ subroutine setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, ver call polympo_getMeshVtxRotLat(mpMesh, nVertices, c_loc(vtxRotLat)) - min = 0 - max = 0 + min = 1000 + max = -1000 do i = 1, nCells sum = 0; do j = 1, nEdgesOnCell(i) @@ -36,7 +36,7 @@ subroutine setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, ver do i = 1, nCells normalizedLat = (elmLat(i) - min) / (max - min) owningProc(i) = normalizedLat * comm_size - owningProc(i) = 0 + if (owningProc(i) == comm_size) owningProc(i) = comm_size-1 end do call polympo_startMeshFill(mpMesh) diff --git a/test/testMigration.cpp b/test/testMigration.cpp index bba49067..736d0030 100644 --- a/test/testMigration.cpp +++ b/test/testMigration.cpp @@ -25,6 +25,7 @@ IntView procWedges(Mesh* mesh, int numElms, int comm_size) { Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){ double normalizedLat = (elmLat(elm) - min(0)) / (max(0) - min(0)); owningProc(elm) = normalizedLat * comm_size; + if (owningProc(elm) == comm_size) owningProc(elm) = comm_size-1; }); return owningProc; } From 6d99d0cc0ab47618c91b01110674d1742020de6b Mon Sep 17 00:00:00 2001 From: Angelyr Date: Tue, 23 Jul 2024 17:32:32 -0700 Subject: [PATCH 24/26] cap normalize --- src/pmpo_wachspressBasis.hpp | 2 +- test/testMigration.cpp | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pmpo_wachspressBasis.hpp b/src/pmpo_wachspressBasis.hpp index a4c64652..77e1e68b 100644 --- a/src/pmpo_wachspressBasis.hpp +++ b/src/pmpo_wachspressBasis.hpp @@ -356,7 +356,7 @@ void sphericalInterpolation(MPMesh& mpMesh){ } }; p_MPs->parallel_for(interpolation, "interpolation"); - pumipic::RecordTime("PolyMPO_Reconstruct", timer.seconds()); + pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds()); } } //namespace polyMPO end diff --git a/test/testMigration.cpp b/test/testMigration.cpp index 736d0030..cd12ad71 100644 --- a/test/testMigration.cpp +++ b/test/testMigration.cpp @@ -23,9 +23,8 @@ IntView procWedges(Mesh* mesh, int numElms, int comm_size) { IntView owningProc("owningProc", numElms); Kokkos::parallel_for("setOwningProc", numElms, KOKKOS_LAMBDA(const int elm){ - double normalizedLat = (elmLat(elm) - min(0)) / (max(0) - min(0)); + double normalizedLat = (elmLat(elm) - min(0)) / (max(0) - min(0)) * .9; owningProc(elm) = normalizedLat * comm_size; - if (owningProc(elm) == comm_size) owningProc(elm) = comm_size-1; }); return owningProc; } @@ -66,7 +65,8 @@ int main(int argc, char* argv[] ) { IntView owningProc = procWedges(mesh, numElms, comm_size); mesh->setMeshEdit(true); mesh->setOwningProc(owningProc); - mpmesh->push(); + for (int i=0; i<5; i++) + mpmesh->push(); auto MPs2Proc = p_MPs.getData(); auto checkPostBackMigrate = PS_LAMBDA(const int& e, const int& mp, const bool& mask) { From 82ba6f4f49864dab3bda7036b31721f6a08a88c5 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 24 Jul 2024 14:11:14 -0700 Subject: [PATCH 25/26] clip normalized --- test/testFortranMPAdvection.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 61ee98db..9478fda1 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -34,9 +34,8 @@ subroutine setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, ver end do do i = 1, nCells - normalizedLat = (elmLat(i) - min) / (max - min) + normalizedLat = (elmLat(i) - min) / (max - min) * .99 owningProc(i) = normalizedLat * comm_size - if (owningProc(i) == comm_size) owningProc(i) = comm_size-1 end do call polympo_startMeshFill(mpMesh) From 50b1893568be2a2011e2652b4e8aa2ab90f39771 Mon Sep 17 00:00:00 2001 From: Angelyr Date: Wed, 14 Aug 2024 12:14:48 -0700 Subject: [PATCH 26/26] fixed migration bug --- src/pmpo_MPMesh.cpp | 6 +++--- test/testFortranMPAdvection.f90 | 33 +++++++++++++++------------------ 2 files changed, 18 insertions(+), 21 deletions(-) diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp index 5922f11d..0e569250 100644 --- a/src/pmpo_MPMesh.cpp +++ b/src/pmpo_MPMesh.cpp @@ -286,14 +286,14 @@ void MPMesh::push(){ bool anyIsMigrating = false; do { CVTTrackingElmCenterBased(); // move to Tgt_XYZ + p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ + p_MPs->updateMPSlice(); // Tgt becomes Cur + p_MPs->updateMPElmID(); //update mpElm IDs slices bool isMigrating = p_MPs->migrate(); anyIsMigrating = getAnyIsMigrating(isMigrating); - p_MPs->updateMPElmID(); //update mpElm IDs slices } while (anyIsMigrating); - p_MPs->updateMPSlice(); // Tgt_XYZ becomes Cur_XYZ - p_MPs->updateMPSlice(); // Tgt becomes Cur pumipic::RecordTime("PolyMPO_push", timer.seconds()); } diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90 index 9478fda1..d0a013f5 100644 --- a/test/testFortranMPAdvection.f90 +++ b/test/testFortranMPAdvection.f90 @@ -2,39 +2,30 @@ module advectionTests contains include "calculateDisplacement.f90" - subroutine setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, verticesOnCell) + subroutine setProcWedges(mpMesh, nCells, comm_size, nEdgesOnCell, verticesOnCell, lonCell) use :: polympo use :: readMPAS use :: iso_c_binding implicit none type(c_ptr) :: mpMesh - integer :: i, j, k, nVertices, nCells, comm_size + integer :: i, j, k, nCells, comm_size integer, dimension(:), pointer :: owningProc, nEdgesOnCell integer, dimension(:,:), pointer :: verticesOnCell - real(kind=MPAS_RKIND), dimension(:), pointer :: elmLat, vtxRotLat - real(kind=MPAS_RKIND) :: normalizedLat, sum, min, max + real(kind=MPAS_RKIND), dimension(:), pointer :: lonCell + real(kind=MPAS_RKIND) :: normalizedLat, min, max - allocate(elmLat(nCells)) allocate(owningProc(nCells)) - allocate(vtxRotLat(nVertices)) - - call polympo_getMeshVtxRotLat(mpMesh, nVertices, c_loc(vtxRotLat)) min = 1000 max = -1000 do i = 1, nCells - sum = 0; - do j = 1, nEdgesOnCell(i) - sum = sum + vtxRotLat(verticesOnCell(j,i)) - end do - elmLat(i) = sum/nEdgesOnCell(i); - if (elmLat(i) < min) min = elmLat(i) - if (elmLat(i) > max) max = elmLat(i) + if (lonCell(i) < min) min = lonCell(i) + if (lonCell(i) > max) max = lonCell(i) end do do i = 1, nCells - normalizedLat = (elmLat(i) - min) / (max - min) * .99 + normalizedLat = (lonCell(i) - min) / (max - min) * .99 owningProc(i) = normalizedLat * comm_size end do @@ -89,7 +80,7 @@ program main real(kind=MPAS_RKIND) :: sphereRadius integer, dimension(:), pointer :: nEdgesOnCell real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex - real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex + real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex, lonCell integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell integer :: numMPs, numMPsCount, numPush integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive @@ -149,6 +140,7 @@ program main print *, "Scale Factor", mpsScaleFactorPerVtx print *, "NUM MPs", numMPs + allocate(lonCell(nCells)) allocate(mpsPerElm(nCells)) allocate(mp2Elm(numMPs)) allocate(isMPActive(numMPs)) @@ -182,6 +174,11 @@ program main yc = yc/nEdgesOnCell(i) zc = zc/nEdgesOnCell(i) + lonCell(i) = atan2(yc,xc) + if (lonCell(i) .le. 0.0_MPAS_RKIND) then ! lon[0,2pi] + lonCell(i) = lonCell(i) + 2.0_MPAS_RKIND*pi + endif + do k = 1, nEdgesOnCell(i) j = verticesOnCell(k,i) @@ -218,7 +215,7 @@ program main call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon)) call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition)) - call setProcWedges(mpMesh, nVertices, nCells, comm_size, nEdgesOnCell, verticesOnCell) + call setProcWedges(mpMesh, nCells, comm_size, nEdgesOnCell, verticesOnCell, lonCell) call runAdvectionTest(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius) call polympo_summarizeTime();