Skip to content

Improve performance of A-take kernels #318

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

Open
julianlitz wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Improve performance of A-take kernels #318

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ef2a35c
Faster Take kernels
julianlitz Aug 16, 2026
96865a0
Add ExtrapolatedSmoother
julianlitz Aug 16, 2026
c875ba7
Format
julianlitz Aug 16, 2026
4e8c545
Fix
julianlitz Aug 16, 2026
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
3 changes: 3 additions & 0 deletions include/Definitions/geometry_helper.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,4 +33,7 @@ KOKKOS_INLINE_FUNCTION void compute_jacobian_elements(const DomainGeometry& doma
/* 1.0 / det(DF) * */
/* [Jtt, -Jrt] */
/* [-Jtr, Jrr] */
/* The stencil coefficient art is scaled by 0.25 to account for the */
/* 9-point stencil used in the finite difference approximation */
art *= 0.25;
}
16 changes: 8 additions & 8 deletions include/DirectSolver/DirectSolverGive/applySymmetryShift.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<do

/* Fill x(i+1,j) */
x(grid.index(i_r + 1, i_theta)) -= -coeff2 * arr * x(grid.index(i_r, i_theta)) /* Left */
+ 0.25 * art * x(grid.index(i_r, i_theta_P1)) /* Top Left */
- 0.25 * art * x(grid.index(i_r, i_theta_M1)); /* Bottom Left */
+ art * x(grid.index(i_r, i_theta_P1)) /* Top Left */
- art * x(grid.index(i_r, i_theta_M1)); /* Bottom Left */

/* --------------------------- */
/* Node next to inner boundary */
Expand All @@ -66,9 +66,9 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<do
/* Fill x(i,j) */
x(grid.index(i_r, i_theta)) -= -coeff1 * arr * x(grid.index(i_r - 1, i_theta)); /* Left */
/* Fill x(i,j-1) */
x(grid.index(i_r, i_theta_M1)) -= +0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Top Left */
x(grid.index(i_r, i_theta_M1)) -= +art * x(grid.index(i_r - 1, i_theta)); /* Top Left */
/* Fill x(i,j+1) */
x(grid.index(i_r, i_theta_P1)) -= -0.25 * art * x(grid.index(i_r - 1, i_theta)); /* Bottom Left */
x(grid.index(i_r, i_theta_P1)) -= -art * x(grid.index(i_r - 1, i_theta)); /* Bottom Left */
}
});
}
Expand Down Expand Up @@ -112,9 +112,9 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<do
/* Fill result(i,j) */
x(grid.index(i_r, i_theta)) -= -coeff2 * arr * x(grid.index(i_r + 1, i_theta)); /* Right */
/* Fill result(i,j-1) */
x(grid.index(i_r, i_theta_M1)) -= -0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Top Right */
x(grid.index(i_r, i_theta_M1)) -= -art * x(grid.index(i_r + 1, i_theta)); /* Top Right */
/* Fill result(i,j+1) */
x(grid.index(i_r, i_theta_P1)) -= +0.25 * art * x(grid.index(i_r + 1, i_theta)); /* Bottom Right */
x(grid.index(i_r, i_theta_P1)) -= +art * x(grid.index(i_r + 1, i_theta)); /* Bottom Right */

/* -------------------------- */
/* Node on the outer boundary */
Expand All @@ -133,8 +133,8 @@ void DirectSolverGive<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<do

/* Fill result(i-1,j) */
x(grid.index(i_r - 1, i_theta)) -= -coeff1 * arr * x(grid.index(i_r, i_theta)) /* Right */
- 0.25 * art * x(grid.index(i_r, i_theta_P1)) /* Top Right */
+ 0.25 * art * x(grid.index(i_r, i_theta_M1)); /* Bottom Right */
- art * x(grid.index(i_r, i_theta_P1)) /* Top Right */
+ art * x(grid.index(i_r, i_theta_M1)); /* Bottom Right */
}
});
}
Expand Down
48 changes: 24 additions & 24 deletions include/DirectSolver/DirectSolverGive/buildSolverMatrix.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,12 +126,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = LeftStencil[StencilPosition::TopRight];
column = top_index;
value = -0.25 * art; /* Top Right */
value = -art; /* Top Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

offset = LeftStencil[StencilPosition::BottomRight];
column = bottom_index;
value = +0.25 * art; /* Bottom Right */
value = +art; /* Bottom Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* Fill matrix row of (i+1,j) */
Expand All @@ -152,12 +152,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = RightStencil[StencilPosition::TopLeft];
column = top_index;
value = +0.25 * art; /* Top Left */
value = +art; /* Top Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

offset = RightStencil[StencilPosition::BottomLeft];
column = bottom_index;
value = -0.25 * art; /* Bottom Left */
value = -art; /* Bottom Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* Fill matrix row of (i,j-1) */
Expand All @@ -178,12 +178,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = BottomStencil[StencilPosition::TopRight];
column = right_index;
value = -0.25 * art; /* Top Right */
value = -art; /* Top Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

offset = BottomStencil[StencilPosition::TopLeft];
column = left_index;
value = +0.25 * art; /* Top Left */
value = +art; /* Top Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* Fill matrix row of (i,j+1) */
Expand All @@ -204,12 +204,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = TopStencil[StencilPosition::BottomRight];
column = right_index;
value = +0.25 * art; /* Bottom Right */
value = +art; /* Bottom Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

offset = TopStencil[StencilPosition::BottomLeft];
column = left_index;
value = -0.25 * art; /* Bottom Left */
value = -art; /* Bottom Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
}
/* -------------------------- */
Expand Down Expand Up @@ -369,12 +369,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = RightStencil[StencilPosition::TopLeft];
column = top_index;
value = +0.25 * art; /* Top Left */
value = +art; /* Top Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

offset = RightStencil[StencilPosition::BottomLeft];
column = bottom_index;
value = -0.25 * art; /* Bottom Left */
value = -art; /* Bottom Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* Fill matrix row of (i,j-1) */
Expand All @@ -395,7 +395,7 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = BottomStencil[StencilPosition::TopRight];
column = right_index;
value = -0.25 * art; /* Top Right */
value = -art; /* Top Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* TopLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
Expand All @@ -418,7 +418,7 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = TopStencil[StencilPosition::BottomRight];
column = right_index;
value = +0.25 * art; /* Bottom Right */
value = +art; /* Bottom Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* BottomLeft REMOVED DUE TO ARTIFICAL 7 POINT STENCIL */
Expand Down Expand Up @@ -512,12 +512,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = LeftStencil[StencilPosition::TopRight];
column = top_index;
value = -0.25 * art; /* Top Right */
value = -art; /* Top Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

offset = LeftStencil[StencilPosition::BottomRight];
column = bottom_index;
value = +0.25 * art; /* Bottom Right */
value = +art; /* Bottom Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
}

Expand All @@ -539,12 +539,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = RightStencil[StencilPosition::TopLeft];
column = top_index;
value = +0.25 * art; /* Top Left */
value = +art; /* Top Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

offset = RightStencil[StencilPosition::BottomLeft];
column = bottom_index;
value = -0.25 * art; /* Bottom Left */
value = -art; /* Bottom Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* Fill matrix row of (i,j-1) */
Expand All @@ -565,14 +565,14 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = BottomStencil[StencilPosition::TopRight];
column = right_index;
value = -0.25 * art; /* Top Right */
value = -art; /* Top Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* REMOVED: Moved to the right hand side to make the matrix symmetric */
if (!DirBC_Interior) {
offset = BottomStencil[StencilPosition::TopLeft];
column = left_index;
value = +0.25 * art; /* Top Left */
value = +art; /* Top Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
}

Expand All @@ -594,14 +594,14 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = TopStencil[StencilPosition::BottomRight];
column = right_index;
value = +0.25 * art; /* Bottom Right */
value = +art; /* Bottom Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* REMOVED: Moved to the right hand side to make the matrix symmetric */
if (!DirBC_Interior) {
offset = TopStencil[StencilPosition::BottomLeft];
column = left_index;
value = -0.25 * art; /* Bottom Left */
value = -art; /* Bottom Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
}
}
Expand Down Expand Up @@ -684,12 +684,12 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = LeftStencil[StencilPosition::TopRight];
column = top_index;
value = -0.25 * art; /* Top Right */
value = -art; /* Top Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

offset = LeftStencil[StencilPosition::BottomRight];
column = bottom_index;
value = 0.25 * art; /* Bottom Right */
value = art; /* Bottom Right */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* Fill matrix row of (i+1,j) */
Expand All @@ -715,7 +715,7 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = BottomStencil[StencilPosition::TopLeft];
column = left_index;
value = 0.25 * art; /* Top Left */
value = art; /* Top Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);

/* Fill matrix row of (i,j+1) */
Expand All @@ -738,7 +738,7 @@ nodeBuildSolverMatrixGive(const int i_r, const int i_theta, const PolarGrid& gri

offset = TopStencil[StencilPosition::BottomLeft];
column = left_index;
value = -0.25 * art; /* Bottom Left */
value = -art; /* Bottom Left */
updateMatrixElement(solver_matrix, ptr, offset, row, column, value);
}
/* ------------------------------------ */
Expand Down
8 changes: 4 additions & 4 deletions include/DirectSolver/DirectSolverTake/applySymmetryShift.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,8 +43,8 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShiftInnerBoundary(Vector<do
const int top = grid.index(i_r, i_theta_P1);

x[center] -= (-coeff1 * (arr[center] + arr[left]) * x[left] /* Left */
- 0.25 * (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+ 0.25 * (art[left] + art[top]) * x[top_left] /* Top Left */
- (art[left] + art[bottom]) * x[bottom_left] /* Bottom Left */
+ (art[left] + art[top]) * x[top_left] /* Top Left */
);
});
}
Expand Down Expand Up @@ -86,8 +86,8 @@ void DirectSolverTake<LevelCacheType>::applySymmetryShiftOuterBoundary(Vector<do
const int top_right = grid.index(i_r + 1, i_theta_P1);

x[center] -= (-coeff2 * (arr[center] + arr[right]) * x[right] /* Right */
+ 0.25 * (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
- 0.25 * (art[right] + art[top]) * x[top_right] /* Top Right */
+ (art[right] + art[bottom]) * x[bottom_right] /* Bottom Right */
- (art[right] + art[top]) * x[top_right] /* Top Right */
);
});
}
Expand Down
24 changes: 12 additions & 12 deletions include/DirectSolver/DirectSolverTake/buildSolverMatrix.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,10 +78,10 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri
- top_value /* Center: (Top) */
);

const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */
const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */
const double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */
const double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */
const double bottom_left_value = -(art(left_index) + art(bottom_index)); /* Bottom Left */
const double bottom_right_value = +(art(right_index) + art(bottom_index)); /* Bottom Right */
const double top_left_value = +(art(left_index) + art(top_index)); /* Top Left */
const double top_right_value = -(art(right_index) + art(top_index)); /* Top Right */

/* Fill matrix row of (i,j) */
row = center_index;
Expand Down Expand Up @@ -204,8 +204,8 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri
- top_value /* Center: (Top) */
);

const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */
const double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */
const double bottom_right_value = +(art(right_index) + art(bottom_index)); /* Bottom Right */
const double top_right_value = -(art(right_index) + art(top_index)); /* Top Right */

/* Fill matrix row of (i,j) */
row = center_index;
Expand Down Expand Up @@ -296,10 +296,10 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri
- top_value /* Center: (Top) */
);

const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */
const double bottom_right_value = +0.25 * (art(right_index) + art(bottom_index)); /* Bottom Right */
const double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */
const double top_right_value = -0.25 * (art(right_index) + art(top_index)); /* Top Right */
const double bottom_left_value = -(art(left_index) + art(bottom_index)); /* Bottom Left */
const double bottom_right_value = +(art(right_index) + art(bottom_index)); /* Bottom Right */
const double top_left_value = +(art(left_index) + art(top_index)); /* Top Left */
const double top_right_value = -(art(right_index) + art(top_index)); /* Top Right */

/* Fill matrix row of (i,j) */
row = center_index;
Expand Down Expand Up @@ -402,8 +402,8 @@ nodeBuildSolverMatrixTake(const int i_r, const int i_theta, const PolarGrid& gri
- top_value /* Center: (Top) */
);

const double bottom_left_value = -0.25 * (art(left_index) + art(bottom_index)); /* Bottom Left */
const double top_left_value = +0.25 * (art(left_index) + art(top_index)); /* Top Left */
const double bottom_left_value = -(art(left_index) + art(bottom_index)); /* Bottom Left */
const double top_left_value = +(art(left_index) + art(top_index)); /* Top Left */

/* Fill matrix row of (i,j) */
row = center_index;
Expand Down
Loading
Loading