Skip to content
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.

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

HPMG: Compute finest level directly on the field #1137

Conversation

AlexanderSinn
Copy link
Member

@AlexanderSinn AlexanderSinn commented Jul 23, 2024

Quite a bit of time in HPMG is spent copying phi between the HPMG array and field array, as the coarsest level has the largest resolution. This PR changes the vcycle so that the finest level is partly directly calculated on the input field array, instead of computing the residual first and then starting at zero. The gsrb_4_residual function is now also used to compute the residual to determine if there should be another vcycle instead of compute_residual, this way gsrb_4_residual can be omitted at the start of the vcycle. This will result in four extra gsrb iterations being done on the finest level, slightly changing CI test results. Additionally, this PR includes some minor cleaning and other smaller HPMG performance improvements.

When using tighter tolerances, the PR and development MG solvers are closer to each other than to the checksum with the lower tolerance (the changes made by this PR are within the tolerance of the MG solver). See
https://github.com/Hi-PACE/hipace/actions/runs/10078089809/job/27862184203?pr=1137
https://github.com/Hi-PACE/hipace/actions/runs/10077289597/job/27859470194?pr=1139

  • Small enough (< few 100s of lines), otherwise it should probably be split into smaller PRs
  • Tested (describe the tests in the PR description)
  • Runs on GPU (basic: the code compiles and run well with the new module)
  • Contains an automated test (checksum and/or comparison with theory)
  • Documented: all elements (classes and their members, functions, namespaces, etc.) are documented
  • Constified (All that can be const is const)
  • Code is clean (no unwanted comments, )
  • Style and code conventions are respected at the bottom of https://github.com/Hi-PACE/hipace
  • Proper label and GitHub project, if applicable

@AlexanderSinn AlexanderSinn added component: fields About 3D fields and slices, field solvers etc. performance optimization, benchmark, profiling, etc. labels Jul 23, 2024
@AlexanderSinn
Copy link
Member Author

On A100 with 4095^2 cells, this provides a 22% performance improvement. @WeiqunZhang could there be accuracy issues from doing gsrb directly on the sol array instead of taking the residual first?

@WeiqunZhang
Copy link
Member

If it converges, I guess it's fine.

@AlexanderSinn
Copy link
Member Author

AlexanderSinn commented Jul 24, 2024

I still need to check if setting the Bx By boundary condition with MR still works the same or if it has to be different now.
... works now.

Copy link
Member

@MaxThevenet MaxThevenet left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this PR! See minor comments below, just for the record.

Comment on lines -728 to +745
if (lev!=0 && (slicemf.box(0).length(0) % 2 == 0)) {
// cell centered MG solve: no ghost cells, put boundary condition into source term
// node centered MG solve: one ghost cell, use boundary condition from there
m_fields.SetBoundaryCondition(m_3D_geom, lev, which_slice, "Bx",
m_fields.getField(lev, which_slice, "Sy"),
amrex::IntVect{0, 0, 0}, 0.5, 8./3.);
m_fields.SetBoundaryCondition(m_3D_geom, lev, which_slice, "By",
m_fields.getField(lev, which_slice, "Sx"),
amrex::IntVect{0, 0, 0}, 0.5, 8./3.);
if (lev!=0) {
if (slicemf.box(0).length(0) % 2 == 0) {
// cell centered MG solve:
m_fields.SetBoundaryCondition(m_3D_geom, lev, which_slice, "Bx",
m_fields.getField(lev, which_slice, "Sy"),
amrex::IntVect{0, 0, 0}, 0.5, 8./3.);
m_fields.SetBoundaryCondition(m_3D_geom, lev, which_slice, "By",
m_fields.getField(lev, which_slice, "Sx"),
amrex::IntVect{0, 0, 0}, 0.5, 8./3.);
} else {
// node centered MG solve:
m_fields.SetBoundaryCondition(m_3D_geom, lev, which_slice, "Bx",
m_fields.getField(lev, which_slice, "Sy"),
amrex::IntVect{0, 0, 0}, 1., 1.);
m_fields.SetBoundaryCondition(m_3D_geom, lev, which_slice, "By",
m_fields.getField(lev, which_slice, "Sx"),
amrex::IntVect{0, 0, 0}, 1., 1.);
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you say a few words on this change?

Comment on lines +823 to +828

// interpolate Bx and By to lev from lev-1 in the ghost cells
m_fields.LevelUpBoundary(m_3D_geom, lev, which_slice, "Bx",
m_fields.m_slices_nguards, m_fields.m_poisson_nguards);
m_fields.LevelUpBoundary(m_3D_geom, lev, which_slice, "By",
m_fields.m_slices_nguards, m_fields.m_poisson_nguards);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also on this one?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because the nodal MG solver doesn't read in the ghost cells anymore, the level >=1 boundary condition has to come from SetBoundaryCondition instead of LevelUpBoundary.

@MaxThevenet MaxThevenet merged commit 008127e into Hi-PACE:development Jul 31, 2024
10 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
component: fields About 3D fields and slices, field solvers etc. performance optimization, benchmark, profiling, etc.
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants