Skip to content

Merge Dn/rotated linear to ac/reconstruction #67

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
onkarsahni merged 17 commits into from
Nov 12, 2024
Merged

Merge Dn/rotated linear to ac/reconstruction #67

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
5f6d963
Option to reverse lat lon increment
dhyan1272 Sep 18, 2024
91b3b9a
Check difference initial and final pos
dhyan1272 Sep 18, 2024
c5792f0
Fixed spherical interp bug, quantiatative push ctest
dhyan1272 Sep 19, 2024
c26a055
Assert statement added
dhyan1272 Sep 20, 2024
f64c2ed
Linear Reconstruction
dhyan1272 Oct 2, 2024
6292949
Linear Construction working, cleaner
dhyan1272 Oct 15, 2024
3dddd47
Cleaned, Linear Reconstruction with regularization
dhyan1272 Nov 4, 2024
333d2de
Some more cleaning
dhyan1272 Nov 4, 2024
644bf31
Added rotated cordinates, added linear reconstruction test
dhyan1272 Nov 7, 2024
726f00f
Merged with ac/reconstruction
dhyan1272 Nov 7, 2024
d90e369
Some more cleaning
dhyan1272 Nov 7, 2024
1da2366
Stricter tolerance
dhyan1272 Nov 7, 2024
e26dc5b
Linear Reconstruction test commented
dhyan1272 Nov 7, 2024
8333764
PR changes 1
dhyan1272 Nov 11, 2024
c7dfb27
Linear rec with MPs as Cell Centres
dhyan1272 Nov 11, 2024
8de3f3f
PR Changes 2
dhyan1272 Nov 11, 2024
dae686b
PR Change 3
dhyan1272 Nov 12, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
89 changes: 80 additions & 9 deletions src/pmpo_MPMesh_assembly.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ void MPMesh::assemblyElm0() {
if(mask) { //if material point is 'active'/'enabled'
Kokkos::atomic_add(&mpsPerElm(elm),1);
for(int j=0;j<numEntries;j++){
Kokkos::atomic_add(&meshField(elm,j), mpData(mp,0));
Kokkos::atomic_add(&meshField(elm,j), mpData(mp,j));
}
}
};
Expand All @@ -100,29 +100,100 @@ void MPMesh::assemblyElm0() {

template <MeshFieldIndex meshFieldIndex>
void MPMesh::assemblyVtx1() {
constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
//Mesh Information
auto elm2VtxConn = p_mesh->getElm2VtxConn();
auto mpData = p_MPs->getData<mpfIndex>();
const int numEntries = mpSliceToNumEntries<mpfIndex>();
int numVtx = p_mesh->getNumVertices();
auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();

//Mesh Field
constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
const int numEntries = mpSliceToNumEntries<mpfIndex>();
p_mesh->fillMeshField<meshFieldIndex>(numVtx, numEntries, 0.0);
auto meshField = p_mesh->getMeshField<meshFieldIndex>();

//Material Points
auto mpData = p_MPs->getData<mpfIndex>();
auto weight = p_MPs->getData<MPF_Basis_Vals>();
auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();

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

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

//Earth Radius
double radius = 1.0;
if(p_mesh->getGeomType() == geom_spherical_surf)
radius=p_mesh->getSphereRadius();

//Assemble matrix for each vertex
auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
if(mask) { //if material point is 'active'/'enabled'
int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element
for(int i=0; i<nVtxE; i++){
int vID = elm2VtxConn(elm,i+1)-1; //vID = vertex id
double fieldComponentVal;
for(int j=0;j<numEntries;j++){
fieldComponentVal = mpData(mp,j);
Kokkos::atomic_add(&meshField(vID,j),fieldComponentVal);
}
double w_vtx=weight(mp,i);

double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius,
(vtxCoords(vID,1) - mpPositions(mp,1))/radius,
(vtxCoords(vID,2) - mpPositions(mp,2))/radius};
for (int k=0; k<vec4d_nEntries; k++)
for (int l=0; l<vec4d_nEntries; l++)
Kokkos::atomic_add(&VtxMatrices(vID,k,l), CoordDiffs[k] * CoordDiffs[l] * w_vtx);
}
}
};
p_MPs->parallel_for(assemble, "assembly");

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

//double f_norm=A.frobeniusNorm();
double A_trace = A.trace();
Matrix4d A_regularized = {v0, v1, v2, v3};
A_regularized.addToDiag(A_trace*1e-8);

double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0};
CholeskySolve4d_UnitRHS(A_regularized, coeff);
for (int i=0; i<vec4d_nEntries; i++)
VtxCoeffs(vtx,i)=coeff[i];
});

//Reconstruct
auto reconstruct = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
if(mask) { //if material point is 'active'/'enabled'
int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element
for(int i=0; i<nVtxE; i++){
int vID = elm2VtxConn(elm,i+1)-1;
double w_vtx=weight(mp,i);
double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius,
(vtxCoords(vID,1) - mpPositions(mp,1))/radius,
(vtxCoords(vID,2) - mpPositions(mp,2))/radius};

auto factor = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] +
VtxCoeffs(vID,2)*CoordDiffs[2] +
VtxCoeffs(vID,3)*CoordDiffs[3]);

for (int k=0; k<numEntries; k++){
auto val = factor*mpData(mp,k);
Kokkos::atomic_add(&reconVals(vID,k), val);
}
}
}
};
p_MPs->parallel_for(reconstruct, "reconstruct");

Kokkos::parallel_for("assigning", numVtx, KOKKOS_LAMBDA(const int vtx){
for(int k=0; k<numEntries; k++)
meshField(vtx, k) = reconVals(vtx,k);
});
}

template <MeshFieldIndex meshFieldIndex>
Expand Down
19 changes: 10 additions & 9 deletions src/pmpo_materialPoints.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -211,25 +211,26 @@ class MaterialPoints {
auto tgtPosXYZ = MPs->get<MPF_Tgt_Pos_XYZ>();
auto rotLatLonIncr = MPs->get<MPF_Rot_Lat_Lon_Incr>();

auto is_rotated = getRotatedFlag();
auto updateRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
if(mask){
auto rotLat = curPosRotLatLon(mp,0) + rotLatLonIncr(mp,0); // phi
auto rotLon = curPosRotLatLon(mp,1) + rotLatLonIncr(mp,1); // lambda
auto rotLon = curPosRotLatLon(mp,1) + rotLatLonIncr(mp,1); // lambda
tgtPosRotLatLon(mp,0) = rotLat;
tgtPosRotLatLon(mp,1) = rotLon;
auto geoLat = rotLat;
auto geoLon = rotLon;
tgtPosRotLatLon(mp,0) = geoLat;
tgtPosRotLatLon(mp,1) = geoLon;
if(is_rotated){
auto xyz_rot = xyz_from_lat_lon(rotLat, rotLon, radius);
auto xyz_geo = grid_rotation_backward(xyz_rot);
lat_lon_from_xyz(geoLat, geoLon, xyz_geo, radius);
}
// x = cosLon cosLat, y = sinLon cosLat, z = sinLat
tgtPosXYZ(mp,0) = radius * std::cos(geoLon) * std::cos(geoLat);
tgtPosXYZ(mp,1) = radius * std::sin(geoLon) * std::cos(geoLat);
tgtPosXYZ(mp,2) = radius * std::sin(geoLat);
tgtPosXYZ(mp,2) = radius * std::sin(geoLat);
}
};
if(isRotatedFlag){
//TODO rotation lat lon calc
fprintf(stderr, "rotational lat lon in MP is not support yet!");
PMT_ALWAYS_ASSERT(false);
}
ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude");
pumipic::RecordTime("PolyMPO_updateRotLatLonAndXYZ2Tgt", timer.seconds());
}
Expand Down
Loading