Skip to content

Commit 8de3f3f

Browse files
dhyan1272dhyan1272
committed
PR Changes 2
1 parent c7dfb27 commit 8de3f3f

File tree

2 files changed

+39
-35
lines changed

2 files changed

+39
-35
lines changed

‎src/pmpo_MPMesh_assembly.hpp‎

Lines changed: 29 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -69,7 +69,6 @@ void MPMesh::assemblyVtx0(){
6969

7070
template <MeshFieldIndex meshFieldIndex>
7171
void MPMesh::assemblyElm0() {
72-
7372
Kokkos::Timer timer;
7473
constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
7574
auto mpData = p_MPs->getData<mpfIndex>();
@@ -84,7 +83,7 @@ void MPMesh::assemblyElm0() {
8483
if(mask) { //if material point is 'active'/'enabled'
8584
Kokkos::atomic_add(&mpsPerElm(elm),1);
8685
for(int j=0;j<numEntries;j++){
87-
Kokkos::atomic_add(&meshField(elm,j), mpData(mp,0));
86+
Kokkos::atomic_add(&meshField(elm,j), mpData(mp,j));
8887
}
8988
}
9089
};
@@ -101,7 +100,6 @@ void MPMesh::assemblyElm0() {
101100

102101
template <MeshFieldIndex meshFieldIndex>
103102
void MPMesh::assemblyVtx1() {
104-
105103
//Mesh Information
106104
auto elm2VtxConn = p_mesh->getElm2VtxConn();
107105
int numVtx = p_mesh->getNumVertices();
@@ -117,16 +115,17 @@ void MPMesh::assemblyVtx1() {
117115
auto mpData = p_MPs->getData<mpfIndex>();
118116
auto weight = p_MPs->getData<MPF_Basis_Vals>();
119117
auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
120-
auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
121118

122119
//Matrix for each vertex
123-
Kokkos::View<double*[4][4]> VtxMatrices("VtxMatrices", p_mesh->getNumVertices());
120+
Kokkos::View<double*[vec4d_nEntries][vec4d_nEntries]> VtxMatrices("VtxMatrices", p_mesh->getNumVertices());
124121

125122
//Reconstructed values
126-
Kokkos::View<double*> reconVals("meshField", p_mesh->getNumVertices());
123+
Kokkos::View<double**> reconVals("meshField", p_mesh->getNumVertices(), numEntries);
127124

128125
//Earth Radius
129-
double radius = p_mesh->getSphereRadius();
126+
double radius = 1.0;
127+
if(p_mesh->getGeomType() == geom_spherical_surf)
128+
radius=p_mesh->getSphereRadius();
130129

131130
//Assemble matrix for each vertex
132131
auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
@@ -136,33 +135,35 @@ void MPMesh::assemblyVtx1() {
136135
int vID = elm2VtxConn(elm,i+1)-1; //vID = vertex id
137136
double w_vtx=weight(mp,i);
138137

139-
double CoordDiffs[4] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius, (vtxCoords(vID,1) - mpPositions(mp,1))/radius,
140-
(vtxCoords(vID,2) - mpPositions(mp,2))/radius}; //add to the matrix
141-
for (int k=0; k<4; k++)
142-
for(int l=0;l<4;l++)
138+
double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius,
139+
(vtxCoords(vID,1) - mpPositions(mp,1))/radius,
140+
(vtxCoords(vID,2) - mpPositions(mp,2))/radius};
141+
for (int k=0; k<vec4d_nEntries; k++)
142+
for (int l=0; l<vec4d_nEntries; l++)
143143
Kokkos::atomic_add(&VtxMatrices(vID,k,l), CoordDiffs[k] * CoordDiffs[l] * w_vtx);
144144
}
145145
}
146146
};
147147
p_MPs->parallel_for(assemble, "assembly");
148148

149149
//Solve Ax=b for each vertex
150-
Kokkos::View<double*[4]> VtxCoeffs("VtxMatrices", p_mesh->getNumVertices());
150+
Kokkos::View<double*[vec4d_nEntries]> VtxCoeffs("VtxMatrices", p_mesh->getNumVertices());
151151
Kokkos::parallel_for("solving Ax=b", numVtx, KOKKOS_LAMBDA(const int vtx){
152152
Vec4d v0 = {VtxMatrices(vtx,0,0), VtxMatrices(vtx,0,1), VtxMatrices(vtx,0,2), VtxMatrices(vtx,0,3)};
153153
Vec4d v1 = {VtxMatrices(vtx,1,0), VtxMatrices(vtx,1,1), VtxMatrices(vtx,1,2), VtxMatrices(vtx,1,3)};
154154
Vec4d v2 = {VtxMatrices(vtx,2,0), VtxMatrices(vtx,2,1), VtxMatrices(vtx,2,2), VtxMatrices(vtx,2,3)};
155155
Vec4d v3 = {VtxMatrices(vtx,3,0), VtxMatrices(vtx,3,1), VtxMatrices(vtx,3,2), VtxMatrices(vtx,3,3)};
156-
Matrix A = {v0,v1,v2,v3};
156+
Matrix4d A = {v0,v1,v2,v3};
157157

158158
//double f_norm=A.frobeniusNorm();
159159
double A_trace = A.trace();
160-
Matrix A_new = {v0, v1, v2, v3};
160+
Matrix4d A_new = {v0, v1, v2, v3};
161161
A_new.regularize(A_trace*1e-8);
162162

163-
double coeff[4]={0.0, 0.0, 0.0, 0.0};
164-
CholeskySolve(A_new, coeff);
165-
for (int i=0; i<4; i++) VtxCoeffs(vtx,i)=coeff[i];
163+
double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0};
164+
CholeskySolve4d(A_new, coeff);
165+
for (int i=0; i<vec4d_nEntries; i++)
166+
VtxCoeffs(vtx,i)=coeff[i];
166167
});
167168

168169
//Reconstruct
@@ -172,18 +173,23 @@ void MPMesh::assemblyVtx1() {
172173
for(int i=0; i<nVtxE; i++){
173174
int vID = elm2VtxConn(elm,i+1)-1;
174175
double w_vtx=weight(mp,i);
175-
double CoordDiffs[4] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius, (vtxCoords(vID,1) - mpPositions(mp,1))/radius,
176-
(vtxCoords(vID,2) - mpPositions(mp,2))/radius};
177-
auto val = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] + VtxCoeffs(vID,2)*CoordDiffs[2] +
178-
VtxCoeffs(vID,3)*CoordDiffs[3])*mpData(mp,0) ;
179-
Kokkos::atomic_add(&reconVals(vID), val);
176+
double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius,
177+
(vtxCoords(vID,1) - mpPositions(mp,1))/radius,
178+
(vtxCoords(vID,2) - mpPositions(mp,2))/radius};
179+
for (int k=0; k<numEntries; k++){
180+
auto val = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] +
181+
VtxCoeffs(vID,2)*CoordDiffs[2] +
182+
VtxCoeffs(vID,3)*CoordDiffs[3])*mpData(mp,k);
183+
Kokkos::atomic_add(&reconVals(vID,k), val);
184+
}
180185
}
181186
}
182187
};
183188
p_MPs->parallel_for(reconstruct, "reconstruct");
184189

185190
Kokkos::parallel_for("assigning", numVtx, KOKKOS_LAMBDA(const int vtx){
186-
meshField(vtx, 0) = reconVals(vtx);
191+
for(int k=0; k<numEntries; k++)
192+
meshField(vtx, k) = reconVals(vtx,k);
187193
});
188194
}
189195

‎src/pmpo_utils.hpp‎

Lines changed: 10 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ typedef double doubleSclr_t[1];
3535
class Vec2d;
3636
class Vec3d;
3737
class Vec4d;
38-
class Matrix;
38+
class Matrix4d;
3939

4040
using Vec2dView = Kokkos::View<Vec2d*>;
4141
using Vec3dView = Kokkos::View<Vec3d*>;
@@ -210,20 +210,18 @@ class Vec4d{
210210
}
211211
};
212212

213-
class Matrix {
213+
class Matrix4d {
214214

215215
private:
216216

217217
double** data_;
218-
int rows_;
219-
int cols_;
218+
int rows_=4;
219+
int cols_=4;
220220

221221
public:
222222

223223
KOKKOS_INLINE_FUNCTION
224-
Matrix(Vec4d v0, Vec4d v1, Vec4d v2, Vec4d v3){
225-
rows_ = 4;
226-
cols_ = 4;
224+
Matrix4d(Vec4d v0, Vec4d v1, Vec4d v2, Vec4d v3){
227225
data_ = new double*[rows_];
228226

229227
for (int i=0; i<rows_; i++){
@@ -240,7 +238,7 @@ class Matrix {
240238

241239
//Destructor
242240
KOKKOS_INLINE_FUNCTION
243-
~Matrix(){
241+
~Matrix4d(){
244242
for (int i=0; i<4; i++){
245243
delete[] data_[i];
246244
}
@@ -253,8 +251,8 @@ class Matrix {
253251

254252
//Scalar division operator overloading
255253
KOKKOS_INLINE_FUNCTION
256-
Matrix operator/(double scalar) const {
257-
Matrix result(Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0));
254+
Matrix4d operator/(double scalar) const {
255+
Matrix4d result(Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0), Vec4d(0, 0, 0, 0));
258256
for (int i = 0; i < rows_; i++) {
259257
for (int j = 0; j < cols_; j++) {
260258
result(i, j) = data_[i][j] / scalar;
@@ -305,7 +303,7 @@ void initArray(Vec2d* arr, int n, Vec2d fill){
305303

306304

307305
KOKKOS_INLINE_FUNCTION
308-
void CholeskySolve(Matrix& A, double* x){
306+
void CholeskySolve4d(Matrix4d& A, double* x){
309307

310308
double a_00=A(0,0);
311309
if (A(0,0)==0){
@@ -370,7 +368,7 @@ void CholeskySolve(Matrix& A, double* x){
370368
}
371369

372370
KOKKOS_INLINE_FUNCTION
373-
void QR_decomp(Matrix& A, Matrix& Q, Matrix &R){
371+
void QR_decomp4d(Matrix4d& A, Matrix4d& Q, Matrix4d& R){
374372
int m_size=4;
375373
for (int i=0; i<m_size; i++){
376374
Vec4d A_column;

0 commit comments

Comments
 (0)