Skip to content

Commit bc2061f

Browse files
dhyan1272dhyan1272
committed
Grid Solve single proc smallest mesh check
1 parent ee2f65e commit bc2061f

6 files changed

Lines changed: 315 additions & 5 deletions

File tree

‎src/pmpo_MPMesh.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -94,6 +94,7 @@ void MPMesh::calculateStress(){
9494

9595
void MPMesh::calculateStressDivergence(){
9696

97+
Kokkos::Timer timer;
9798
int self, numProcsTot;
9899
MPI_Comm comm = p_MPs->getMPIComm();
99100
MPI_Comm_rank(comm, &self);
@@ -183,6 +184,8 @@ void MPMesh::calculateStressDivergence(){
183184
stressDivergence(vtx, 0) = ramp * stress_divUV(vtx, 0) + (1 - ramp) * divUV_edge(vtx, 0) * invM ;
184185
stressDivergence(vtx, 1) = ramp * stress_divUV(vtx, 1) + (1 - ramp) * divUV_edge(vtx, 1) * invM ;
185186
});
187+
188+
pumipic::RecordTime("Stress_Divergence_Reconstruction" + std::to_string(self), timer.seconds());
186189
}
187190

188191
void MPMesh::calcBasis() {

‎src/pmpo_c.cpp‎

Lines changed: 119 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1371,11 +1371,26 @@ void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* ar
13711371
Kokkos::deep_copy(solveStress, h_solveStress);
13721372
}
13731373

1374+
void polympo_setSolveVelocityMesh_f(MPMesh_ptr p_mpmesh, const int nVertices, int* array){
1375+
//chech validity
1376+
checkMPMeshValid(p_mpmesh);
1377+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1378+
1379+
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
1380+
//copy the host array to the device
1381+
auto solveVelocity = p_mesh->getMeshField<polyMPO::MeshF_SolveVelocity>();
1382+
auto h_solveVelocity = Kokkos::create_mirror_view(solveVelocity);
1383+
for(int i=0; i<nVertices; i++)
1384+
h_solveVelocity(i) = array[i];
1385+
Kokkos::deep_copy(solveVelocity, h_solveVelocity);
1386+
}
1387+
13741388
void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh){
13751389
//chech validity
13761390
checkMPMeshValid(p_mpmesh);
13771391
((polyMPO::MPMesh*)p_mpmesh) -> calculateStressDivergence();
1378-
}
1392+
1393+
}
13791394

13801395
void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray){
13811396
//chech validity
@@ -1394,6 +1409,109 @@ void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, dou
13941409
}
13951410
}
13961411

1412+
void polympo_setTotalMassVtx_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){
1413+
checkMPMeshValid(p_mpmesh);
1414+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1415+
1416+
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
1417+
//copy the host array to the device
1418+
auto totalMassVtx = p_mesh->getMeshField<polyMPO::MeshF_TotalMassVtx>();
1419+
auto h_totalMassVtx = Kokkos::create_mirror_view(totalMassVtx);
1420+
for(int i=0; i<nVertices; i++)
1421+
h_totalMassVtx(i, 0) = array[i];
1422+
Kokkos::deep_copy(totalMassVtx, h_totalMassVtx);
1423+
}
1424+
1425+
void polympo_set_airStress_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray){
1426+
//check mpMesh is valid
1427+
checkMPMeshValid(p_mpmesh);
1428+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1429+
1430+
//check the size
1431+
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
1432+
1433+
//copy the host array to the device
1434+
auto airStress = p_mesh->getMeshField<polyMPO::MeshF_AirStress>();
1435+
auto h_airStress = Kokkos::create_mirror_view(airStress);
1436+
for(int i=0; i<nVertices; i++){
1437+
h_airStress(i,0) = uArray[i];
1438+
h_airStress(i,1) = vArray[i];
1439+
}
1440+
Kokkos::deep_copy(airStress, h_airStress);
1441+
}
1442+
1443+
void polympo_set_surfaceTiltForce_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray){
1444+
//check mpMesh is valid
1445+
checkMPMeshValid(p_mpmesh);
1446+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1447+
1448+
//check the size
1449+
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
1450+
1451+
//copy the host array to the device
1452+
auto surfTiltForce = p_mesh->getMeshField<polyMPO::MeshF_SurfaceTilt>();
1453+
auto h_surfTiltForce = Kokkos::create_mirror_view(surfTiltForce);
1454+
for(int i=0; i<nVertices; i++){
1455+
h_surfTiltForce(i,0) = uArray[i];
1456+
h_surfTiltForce(i,1) = vArray[i];
1457+
}
1458+
Kokkos::deep_copy(surfTiltForce, h_surfTiltForce);
1459+
}
1460+
1461+
void polympo_set_totalMassVertexfVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){
1462+
//check mpMesh is valid
1463+
checkMPMeshValid(p_mpmesh);
1464+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1465+
1466+
//check the size
1467+
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
1468+
1469+
//copy the host array to the device
1470+
auto totalMassFVtx = p_mesh->getMeshField<polyMPO::MeshF_TotalMassFVtx>();
1471+
auto h_totalMassFVtx = Kokkos::create_mirror_view(totalMassFVtx);
1472+
for(int i=0; i<nVertices; i++)
1473+
h_totalMassFVtx(i,0) = array[i];
1474+
Kokkos::deep_copy(totalMassFVtx, h_totalMassFVtx);
1475+
}
1476+
1477+
void polympo_set_oceanStress_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray){
1478+
//check mpMesh is valid
1479+
checkMPMeshValid(p_mpmesh);
1480+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1481+
//check the size
1482+
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
1483+
1484+
//copy the host array to the device
1485+
auto oceanStress = p_mesh->getMeshField<polyMPO::MeshF_OceanStress>();
1486+
auto h_oceanStress = Kokkos::create_mirror_view(oceanStress);
1487+
for(int i=0; i<nVertices; i++){
1488+
h_oceanStress(i,0) = uArray[i];
1489+
h_oceanStress(i,1) = vArray[i];
1490+
}
1491+
Kokkos::deep_copy(oceanStress, h_oceanStress);
1492+
}
1493+
1494+
void polympo_set_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array){
1495+
//check mpMesh is valid
1496+
checkMPMeshValid(p_mpmesh);
1497+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1498+
//check the size
1499+
PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
1500+
1501+
//copy the host array to the device
1502+
auto oceanStressCoeff = p_mesh->getMeshField<polyMPO::MeshF_OceanStressCoeff>();
1503+
auto h_oceanStressCoeff = Kokkos::create_mirror_view(oceanStressCoeff);
1504+
for(int i=0; i<nVertices; i++)
1505+
h_oceanStressCoeff(i,0) = array[i];
1506+
1507+
Kokkos::deep_copy(oceanStressCoeff, h_oceanStressCoeff);
1508+
}
1509+
1510+
void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh){
1511+
//Temporary grid Solve after calculateDivergence
1512+
auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
1513+
p_mesh->gridSolveGPU();
1514+
}
13971515

13981516
//Advection Calcualtions
13991517
void polympo_push_f(MPMesh_ptr p_mpmesh){

‎src/pmpo_c.h‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -102,9 +102,16 @@ void polyMPO_setTanLatVertexRotatedOverRadius_f(MPMesh_ptr p_mpmesh, const int n
102102
void polympo_setElasticTimeStep_f(MPMesh_ptr p_mpmesh, const double elasticTimeStep);
103103
void polympo_setDynamicTimeStep_f(MPMesh_ptr p_mpmesh, const double dynamicTimeStep);
104104
void polympo_setSolveStressMesh_f(MPMesh_ptr p_mpmesh, const int nCells, int* array);
105+
void polympo_setSolveVelocityMesh_f(MPMesh_ptr p_mpmesh, const int nVertices, int* array);
105106
void polympo_calculateStressDivergence_f(MPMesh_ptr p_mpmesh);
106107
void polympo_getStressDivergence_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray);
107-
108+
void polympo_setTotalMassVtx_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
109+
void polympo_set_airStress_f(MPMesh_ptr mpmesh, const int nVertices, double* uArray,double* vArray);
110+
void polympo_set_surfaceTiltForce_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray);
111+
void polympo_set_totalMassVertexfVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
112+
void polympo_set_oceanStress_f(MPMesh_ptr p_mpmesh, const int nVertices, double* uArray, double* vArray);
113+
void polympo_set_oceanStressCoefficient_f(MPMesh_ptr p_mpmesh, const int nVertices, double* array);
114+
void polympo_velocity_grid_solve_f(MPMesh_ptr p_mpmesh);
108115

109116
// Advection calculations
110117
void polympo_push_f(MPMesh_ptr p_mpmesh);

‎src/pmpo_fortran.f90‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -976,6 +976,14 @@ subroutine polympo_setSolveStressMesh(mpMesh, nCells, array) &
976976
type(c_ptr), value :: array
977977
end subroutine
978978

979+
subroutine polympo_setSolveVelocityMesh(mpMesh, nVertices, array) &
980+
bind(C, NAME='polympo_setSolveVelocityMesh_f')
981+
use :: iso_c_binding
982+
type(c_ptr), value :: mpMesh
983+
integer(c_int), value :: nVertices
984+
type(c_ptr), value :: array
985+
end subroutine
986+
979987
subroutine polympo_calculateStressDivergence(mpMesh) &
980988
bind(C, NAME='polympo_calculateStressDivergence_f')
981989
use :: iso_c_binding
@@ -990,6 +998,61 @@ subroutine polympo_getStressDivergence(mpMesh, nVertices, uArray, vArray) &
990998
type(c_ptr), value :: uArray, vArray
991999
end subroutine
9921000

1001+
subroutine polympo_setTotalMassVtx(mpMesh, nVertices, array) &
1002+
bind(C, NAME='polympo_setTotalMassVtx_f')
1003+
use :: iso_c_binding
1004+
type(c_ptr), value :: mpMesh
1005+
integer(c_int), value :: nVertices
1006+
type(c_ptr), value :: array
1007+
end subroutine
1008+
1009+
subroutine polympo_set_airStress(mpMesh, nVertices, uArray, vArray) &
1010+
bind(C, NAME='polympo_set_airStress_f')
1011+
use :: iso_c_binding
1012+
type(c_ptr), value :: mpMesh
1013+
integer(c_int), value :: nVertices
1014+
type(c_ptr), value :: uArray, vArray
1015+
end subroutine
1016+
1017+
subroutine polympo_set_surfaceTiltForce(mpMesh, nVertices, uArray, vArray) &
1018+
bind(C, NAME='polympo_set_surfaceTiltForce_f')
1019+
use :: iso_c_binding
1020+
type(c_ptr), value :: mpMesh
1021+
integer(c_int), value :: nVertices
1022+
type(c_ptr), value :: uArray, vArray
1023+
end subroutine
1024+
1025+
subroutine polympo_set_totalMassVertexfVertex(mpMesh, nVertices, array) &
1026+
bind(C, NAME='polympo_set_totalMassVertexfVertex_f')
1027+
use :: iso_c_binding
1028+
type(c_ptr), value :: mpMesh
1029+
integer(c_int), value :: nVertices
1030+
type(c_ptr), value :: array
1031+
end subroutine
1032+
1033+
subroutine polympo_set_oceanStress(mpMesh, nVertices, uArray, vArray) &
1034+
bind(C, NAME='polympo_set_oceanStress_f')
1035+
use :: iso_c_binding
1036+
type(c_ptr), value :: mpMesh
1037+
integer(c_int), value :: nVertices
1038+
type(c_ptr), value :: uArray, vArray
1039+
end subroutine
1040+
1041+
subroutine polympo_set_oceanStressCoefficient(mpMesh, nVertices, array) &
1042+
bind(C, NAME='polympo_set_oceanStressCoefficient_f')
1043+
use :: iso_c_binding
1044+
type(c_ptr), value :: mpMesh
1045+
integer(c_int), value :: nVertices
1046+
type(c_ptr), value :: array
1047+
end subroutine
1048+
1049+
subroutine polympo_velocity_grid_solve(mpMesh) &
1050+
bind(C, NAME='polympo_velocity_grid_solve_f')
1051+
use :: iso_c_binding
1052+
type(c_ptr), value :: mpMesh
1053+
end subroutine
1054+
1055+
9931056

9941057
!---------------------------------------------------------------------------
9951058
!> @brief calculate the MPs from given mesh vertices rotational latitude

‎src/pmpo_mesh.cpp‎

Lines changed: 67 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -47,14 +47,33 @@ namespace polyMPO{
4747
auto tanLatVertexRotOverRadiusEntry = meshFields2TypeAndString.at(MeshF_TanLatVertexRotatedOverRadius);
4848
PMT_ALWAYS_ASSERT(tanLatVertexRotOverRadiusEntry.first == MeshFType_VtxBased);
4949
tanLatVertexRotatedOverRadius_ = MeshFView<MeshF_TanLatVertexRotatedOverRadius>(tanLatVertexRotOverRadiusEntry.second, numVtxs_);
50-
50+
5151
auto interiorVertexEntry = meshFields2TypeAndString.at(MeshF_InteriorVertex);
5252
PMT_ALWAYS_ASSERT(interiorVertexEntry.first == MeshFType_VtxBased);
5353
interiorVertex_ = MeshFView<MeshF_InteriorVertex>(interiorVertexEntry.second, numVtxs_);
5454

5555
auto stressDivergenceEntry = meshFields2TypeAndString.at(MeshF_StressDivergence);
5656
PMT_ALWAYS_ASSERT(stressDivergenceEntry.first == MeshFType_VtxBased);
5757
stressDivergence_ = MeshFView<MeshF_StressDivergence>(stressDivergenceEntry.second, numVtxs_);
58+
59+
auto solveVelocityEntry = meshFields2TypeAndString.at(MeshF_SolveVelocity);
60+
PMT_ALWAYS_ASSERT(solveVelocityEntry.first == MeshFType_VtxBased);
61+
solveVelocity_ = MeshFView<MeshF_SolveVelocity>(solveVelocityEntry.second, numVtxs_);
62+
63+
auto totalMassVtxEntry = meshFields2TypeAndString.at(MeshF_TotalMassVtx);
64+
PMT_ALWAYS_ASSERT(totalMassVtxEntry.first == MeshFType_VtxBased);
65+
totalMassVtx_ = MeshFView<MeshF_TotalMassVtx>(totalMassVtxEntry.second, numVtxs_);
66+
67+
//For grid solve, the source terms
68+
airStress_ = MeshFView<MeshF_AirStress>(meshFields2TypeAndString.at(MeshF_AirStress).second, numVtxs_);
69+
70+
surfaceTilt_ = MeshFView<MeshF_SurfaceTilt>(meshFields2TypeAndString.at(MeshF_SurfaceTilt).second, numVtxs_);
71+
72+
totalMassFVtx_ = MeshFView<MeshF_TotalMassFVtx>(meshFields2TypeAndString.at(MeshF_TotalMassFVtx).second, numVtxs_);
73+
74+
oceanStress_ = MeshFView<MeshF_OceanStress>(meshFields2TypeAndString.at(MeshF_OceanStress).second, numVtxs_);
75+
76+
oceanStressCoeff_ = MeshFView<MeshF_OceanStressCoeff>(meshFields2TypeAndString.at(MeshF_OceanStressCoeff).second, numVtxs_);
5877
}
5978

6079
void Mesh::setMeshElmBasedFieldSize(){
@@ -134,4 +153,51 @@ namespace polyMPO{
134153
});
135154
}
136155

156+
void Mesh::gridSolveGPU(){
157+
std::cout<<__FUNCTION__<<std::endl;
158+
//Mesh Fields
159+
int numVertices = getNumVertices();
160+
auto totalMassVtx = getMeshField<MeshF_TotalMassVtx>();
161+
auto totalMassFVtx = getMeshField<MeshF_TotalMassFVtx>();
162+
auto airStress = getMeshField<MeshF_AirStress>();
163+
auto surfaceTiltForce = getMeshField<MeshF_SurfaceTilt>();
164+
auto oceanStress = getMeshField<MeshF_OceanStress>();
165+
auto oceanStressCoeff = getMeshField<MeshF_OceanStressCoeff>();
166+
auto solve_velocity = getMeshField<MeshF_SolveVelocity>();
167+
auto stressDivergence = getMeshField<MeshF_StressDivergence>();
168+
auto velocity = getMeshField<MeshF_Vel>();
169+
auto elasticTimeStep = getElasticTimeStep();
170+
171+
double sinOceanTurningAngle=0.0;
172+
double cosOceanTurningAngle=1.0;
173+
174+
Kokkos::parallel_for("SolveGridVelocity", numVertices, KOKKOS_LAMBDA(const int vtx){
175+
if(solve_velocity(vtx) == 0) return;
176+
double a, b, c, d, s, rhs_u, rhs_v, denom ;
177+
178+
s = (totalMassFVtx(vtx, 0) >= 0) ? 1.0 : -1.0;
179+
a = totalMassVtx(vtx, 0)/elasticTimeStep + oceanStressCoeff(vtx, 0) * cosOceanTurningAngle;
180+
b = -totalMassFVtx(vtx, 0) - oceanStressCoeff(vtx, 0) * sinOceanTurningAngle * s * totalMassFVtx(vtx, 0);
181+
c = totalMassFVtx(vtx, 0) + oceanStressCoeff(vtx, 0) * sinOceanTurningAngle * s * totalMassFVtx(vtx, 0);
182+
d = totalMassVtx(vtx, 0)/elasticTimeStep + oceanStressCoeff(vtx, 0) * cosOceanTurningAngle;
183+
rhs_u = stressDivergence(vtx, 0) + airStress(vtx, 0) + surfaceTiltForce(vtx, 0) + oceanStressCoeff(vtx, 0)*oceanStress(vtx, 0) +
184+
(totalMassVtx(vtx, 0) * velocity(vtx, 0))/elasticTimeStep;
185+
rhs_v = stressDivergence(vtx, 1) + airStress(vtx, 1) + surfaceTiltForce(vtx, 1) + oceanStressCoeff(vtx, 0)*oceanStress(vtx, 1) +
186+
(totalMassVtx(vtx, 0) * velocity(vtx, 1))/elasticTimeStep;
187+
denom = a*d - b*c;
188+
189+
velocity(vtx, 0) = (d*rhs_u-b*rhs_v)/denom;
190+
velocity(vtx, 1) = (a*rhs_v-c*rhs_u)/denom;
191+
192+
//Debug
193+
if(vtx == 10 || vtx == 11){
194+
/*
195+
printf("Vtx %d F: %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e %.15e \n", vtx, totalMassVtx(vtx,0), totalMassFVtx(vtx,0),
196+
airStress(vtx,0), airStress(vtx,1), surfaceTiltForce(vtx,0), surfaceTiltForce(vtx,1), oceanStress(vtx,0), oceanStress(vtx,1),
197+
oceanStressCoeff(vtx, 0));*/
198+
printf("Vtx %d V: %.15e %.15e \n", vtx, velocity(vtx,0), velocity(vtx,1));
199+
}
200+
});
201+
}
202+
137203
} // namespace polyMPO

0 commit comments

Comments
 (0)