@@ -52,10 +52,10 @@ void MPMesh::calculateStrain(){
5252 v22 = v22 + MPsBasisGrads (mp, i*2 + 1 ) * velField (iVertex, 1 );
5353 uTanOverR = uTanOverR + MPsBasis (mp, i) * tanLatVertexRotatedOverRadius (iVertex, 0 ) * velField (iVertex, 0 );
5454 vTanOverR = vTanOverR + MPsBasis (mp, i) * tanLatVertexRotatedOverRadius (iVertex, 0 ) * velField (iVertex, 1 );
55- // Debugging
56- if (MPsAppID (mp)==0 ){
57- // printf("Strain Calc: iVertex %d vel field %.15e %.15e \n", iVertex, velField(iVertex, 0), velField(iVertex, 1));
58- }
55+ /*
56+ if(MPsAppID(mp)==0)
57+ printf("Strain Calc: iVertex %d vel field %.15e %.15e \n", iVertex, velField(iVertex, 0), velField(iVertex, 1));
58+ */
5959 }
6060 MPsStrainRate (mp, 0 ) = v11 - vTanOverR;
6161 MPsStrainRate (mp, 1 ) = v22;
@@ -84,21 +84,28 @@ void MPMesh::calculateStress(){
8484
8585 Vec3d strain_rate (MPsStrainRate (mp, 0 ), MPsStrainRate (mp, 1 ), MPsStrainRate (mp, 2 ));
8686 Vec3d stress (MPsStress (mp, 0 ), MPsStress (mp, 1 ), MPsStress (mp, 2 ));
87- constitutive_evp (strain_rate, stress, MPsIcePressure (mp, 0 ), MPsRepPressure (mp, 0 ), MPsArea (mp, 0 ), elasticTimeStep, dampingTimescale);
87+ // constitutive_evp(strain_rate, stress, MPsIcePressure(mp, 0), MPsRepPressure(mp, 0), MPsArea(mp, 0), elasticTimeStep, dampingTimescale);
88+ constitutive_linear (strain_rate, stress);
8889 for (int m=0 ; m<3 ; m++)
8990 MPsStress (mp, m) = stress[m];
90- // Debugging
91+ /*
9192 if(MPsAppID(mp)==0){
92- // printf("Strain in GPU: %.15e %.15e %.15e\n", MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
93- // printf("Stress in GPU: %.15e %.15e %.15e\n", MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
93+ printf("Strain in GPU: %.15e %.15e %.15e\n", MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
94+ printf("Stress in GPU: %.15e %.15e %.15e\n", MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
9495 }
96+ */
9597 }
9698 };
9799 p_MPs->parallel_for (setMPStress, " setMPStress" );
98100}
99101
100102void MPMesh::calculateStressDivergence (){
101-
103+
104+ int self, numProcsTot;
105+ MPI_Comm comm = p_MPs->getMPIComm ();
106+ MPI_Comm_rank (comm, &self);
107+ MPI_Comm_size (comm, &numProcsTot);
108+
102109 // Mesh Information
103110 auto elm2VtxConn = p_mesh->getElm2VtxConn ();
104111 int numVtx = p_mesh->getNumVertices ();
@@ -125,11 +132,11 @@ void MPMesh::calculateStressDivergence(){
125132 radius=p_mesh->getSphereRadius ();
126133
127134 // Reconstructed the stress
128- Kokkos::View<double *> stress_divU (" stress_divu" , p_mesh->getNumVertices ());
129- Kokkos::View<double *> stress_divV (" stress_divv" , p_mesh->getNumVertices ());
135+ Kokkos::View<doubleSclr_t *> stress_divU (" stress_divu" , p_mesh->getNumVertices ());
136+ Kokkos::View<doubleSclr_t *> stress_divV (" stress_divv" , p_mesh->getNumVertices ());
130137
131- Kokkos::View<double *> divU_edge (" divUedge" , p_mesh->getNumVertices ());
132- Kokkos::View<double *> divV_edge (" divVedge" , p_mesh->getNumVertices ());
138+ Kokkos::View<doubleSclr_t *> divU_edge (" divUedge" , p_mesh->getNumVertices ());
139+ Kokkos::View<doubleSclr_t *> divV_edge (" divVedge" , p_mesh->getNumVertices ());
133140
134141 // Assemble fields for Stress Divergence
135142 auto stress_div = PS_LAMBDA (const int & elm, const int & mp, const int & mask) {
@@ -153,45 +160,53 @@ void MPMesh::calculateStressDivergence(){
153160 VtxCoeffs_new (vID,2 , 2 )*CoordDiffs[2 ] +
154161 VtxCoeffs_new (vID,2 , 3 )*CoordDiffs[3 ]);
155162
156- Kokkos::atomic_add (&stress_divU (vID), factor1 * MPsStress (mp, 0 ) + factor2 * MPsStress (mp, 2 ) -
163+ Kokkos::atomic_add (&stress_divU (vID, 0 ), factor1 * MPsStress (mp, 0 ) + factor2 * MPsStress (mp, 2 ) -
157164 2 * tanLatVertexRotatedOverRadius (vID, 0 ) * factor * MPsStress (mp, 2 ));
158- Kokkos::atomic_add (&stress_divV (vID), factor2 * MPsStress (mp, 1 ) + factor1*MPsStress (mp, 2 ) +
165+ Kokkos::atomic_add (&stress_divV (vID, 0 ), factor2 * MPsStress (mp, 1 ) + factor1*MPsStress (mp, 2 ) +
159166 factor * tanLatVertexRotatedOverRadius (vID, 0 ) * (MPsStress (mp, 0 )-MPsStress (mp, 1 )));
160167
161168
162- Kokkos::atomic_add (&divU_edge (vID), - weight_grads (mp, i*2 + 0 ) * MPsStress (mp, 0 ) - weight_grads (mp, i*2 + 1 ) * MPsStress (mp, 2 ) -
169+ Kokkos::atomic_add (&divU_edge (vID, 0 ), - weight_grads (mp, i*2 + 0 ) * MPsStress (mp, 0 ) - weight_grads (mp, i*2 + 1 ) * MPsStress (mp, 2 ) -
163170 2 * tanLatVertexRotatedOverRadius (vID, 0 ) * w_vtx * MPsStress (mp, 2 ));
164- Kokkos::atomic_add (&divV_edge (vID), - weight_grads (mp, i*2 + 1 ) * MPsStress (mp, 1 ) - weight_grads (mp, i*2 + 0 ) * MPsStress (mp, 2 ) +
171+ Kokkos::atomic_add (&divV_edge (vID, 0 ), - weight_grads (mp, i*2 + 1 ) * MPsStress (mp, 1 ) - weight_grads (mp, i*2 + 0 ) * MPsStress (mp, 2 ) +
165172 w_vtx * tanLatVertexRotatedOverRadius (vID, 0 ) * (MPsStress (mp, 0 )- MPsStress (mp, 1 )));
166173
167174 }
168175 }
169176 };
170- p_MPs->parallel_for (stress_div, " assembly" );
171-
172- // TODO COMMUNICATE THE VERTEX FIELDS
173-
177+ p_MPs->parallel_for (stress_div, " stress_div_assembly" );
178+
179+ // TO DO make to 2 fields instrad of 4 to reduce latency
180+ if (numProcsTot>1 ){
181+ // Takes contribution of halo vertices and adds it in halo cells
182+ communicate_and_take_halo_contributions (stress_divU, numVertices, 1 , 0 , 0 );
183+ communicate_and_take_halo_contributions (stress_divV, numVertices, 1 , 0 , 0 );
184+ communicate_and_take_halo_contributions (divU_edge, numVertices, 1 , 0 , 0 );
185+ communicate_and_take_halo_contributions (divV_edge, numVertices, 1 , 0 , 0 );
186+ // Does it need to transfer the correct values at owned vertices to halo vertices
187+ communicate_and_take_halo_contributions (stress_divU, numVertices, 1 , 1 , 1 );
188+ communicate_and_take_halo_contributions (stress_divV, numVertices, 1 , 1 , 1 );
189+ communicate_and_take_halo_contributions (divU_edge, numVertices, 1 , 1 , 1 );
190+ communicate_and_take_halo_contributions (divV_edge, numVertices, 1 , 1 , 1 );
191+ }
192+
174193 auto stressDivergence = p_mesh->getMeshField <MeshF_StressDivergence>();
175-
176194 Kokkos::parallel_for (" calculate_divergence" , numVtx, KOKKOS_LAMBDA (const int vtx){
177195
178196 double ramp = nearAnEdge_l (vtx);
179197 double invM = 1.0 /vtxMatrixMass_l (vtx);
180198 invM = vtxMatrixMass_l (vtx) >1e-4 ? invM : 0 ;
181199
182- stressDivergence (vtx, 0 ) = ramp * stress_divU (vtx) + (1 - ramp) * divU_edge (vtx) * invM ;
183- stressDivergence (vtx, 1 ) = ramp * stress_divV (vtx) + (1 - ramp) * divV_edge (vtx) * invM ;
184- // Debugging
185- if (vtx == 10 || vtx == 11 || vtx == 1476 || vtx == 1481 ) {
186- // printf("Vtx %d Divergence %.15e %.15e \n", vtx, stressDivergence(vtx, 0), stressDivergence(vtx, 1));
200+ stressDivergence (vtx, 0 ) = ramp * stress_divU (vtx, 0 ) + (1 - ramp) * divU_edge (vtx, 0 ) * invM ;
201+ stressDivergence (vtx, 1 ) = ramp * stress_divV (vtx, 0 ) + (1 - ramp) * divV_edge (vtx, 0 ) * invM ;
202+ /*
203+ if (ent2global( vtx) == 1714 ) {
204+ printf("Vtx %d Divergence %.15e %.15e %.15e %.15e \n", vtx, stressDivergence(vtx, 0), stressDivergence(vtx, 1));
187205 }
206+ */
188207 });
189-
190- // TODO COMMUNICATE THE VERTEX FIELDS
191-
192208}
193209
194-
195210void MPMesh::calcBasis () {
196211 assert (p_mesh->getGeomType () == geom_spherical_surf);
197212
0 commit comments