Skip to content

Conversation

@AlexanderSinn
Copy link
Member

@AlexanderSinn AlexanderSinn commented Feb 9, 2022

based on #723, #724 and #725

To enable Open Boundaries use:

geometry.is_periodic = false false false
fields.extended_solve = true
fields.open_boundary = true

The Poisson equations are solved with open boundaries rather than zero boundaries by integrating a Taylor expansion of the Greens function of the transverse laplacien with the source and using that as boundary conditions of the FFT Poisson solver. Note that this does not apply to the Helmholtz multigrid solver.

Old: development
New: Open Boundary

image
image
image
image

Performance:
amr.n_cell = 1019 1019 2000 and plasma.ppc = 2 2
Open Boundaries are calculated in Fields::SetBoundaryCondition()

--------------------------------------------------------------------------------------------------
Name                                               NCalls  Excl. Min  Excl. Avg  Excl. Max   Max %
--------------------------------------------------------------------------------------------------
DepositCurrent_PlasmaParticleContainer()            25674      33.66      33.66      33.66  24.92%
UpdateForcePushParticles_Plasma(tp  )               23673      24.52      24.52      24.52  18.16%
Fields::SetBoundaryCondition()                      53346      15.54      15.54      15.54  11.51%
UpdateForcePushParticles_Plasma(  u )               23673      14.26      14.26      14.26  10.56%
Hipace::PredictorCorrectorLoopToSolveBxBy()          2000      5.857      5.857      5.857   4.34%
AnyDST::C2Rfft()                                   213384      4.769      4.769      4.769   3.53%
AnyDST::ToSine()                                   213384      4.409      4.409      4.409   3.27%
AnyDST::Transpose()                                213384      3.606      3.606      3.606   2.67%
AnyDST::ToComplex()                                213384      3.414      3.414      3.414   2.53%
Fields::ComputeRelBFieldError()                     25673      2.906      2.906      2.906   2.15%
MultiFab::LinComb()                                 98692      2.676      2.676      2.676   1.98%
FFTPoissonSolverDirichlet::SolvePoissonEquation()   53346      2.419      2.419      2.419   1.79%

amr.n_cell = 2043 2043 3000 and plasma.ppc = 1 1:

--------------------------------------------------------------------------------------------------
Name                                               NCalls  Excl. Min  Excl. Avg  Excl. Max   Max %
--------------------------------------------------------------------------------------------------
Fields::SetBoundaryCondition()                      92368      80.64      80.64      80.64  18.90%
DepositCurrent_PlasmaParticleContainer()            44685       66.7       66.7       66.7  15.63%
UpdateForcePushParticles_Plasma(  u )               41684      47.39      47.39      47.39  11.11%
UpdateForcePushParticles_Plasma(tp  )               41684      43.24      43.24      43.24  10.13%
AnyDST::C2Rfft()                                   369472      23.63      23.63      23.63   5.54%
AnyDST::ToSine()                                   369472      23.09      23.09      23.09   5.41%
AnyDST::Transpose()                                369472         23         23         23   5.39%
AnyDST::ToComplex()                                369472       19.6       19.6       19.6   4.59%
MultiFab::LinComb()                                172736      14.69      14.69      14.69   3.44%
FFTPoissonSolverDirichlet::SolvePoissonEquation()   92368      13.22      13.22      13.22   3.10%

@SeverinDiederichs
Copy link
Member

SeverinDiederichs commented Mar 1, 2022

image

This is an offset with (x,y) = (2,2).

Looks great!

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 inlined. Besides,

  • Could you see with @SeverinDiederichs to add the offset beam in vacuum test in CI?
  • Could you share profiling data for 1 ppc and a transverse box size of 2048^2?

// ignore everything outside of 95% the min radius as the Taylor expansion only converges
// outside of a circular patch containing the sources, i.e. the sources can't be further
// from the center than the closest boundary as it would be the case in the corners
const amrex::Real cutoff_sq = pow<2>(0.95_rt * radius * scale);
Copy link
Member

Choose a reason for hiding this comment

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

Isn't 95% arbitrary? Should it be e.g. min radius minus 3 cells or so?

Copy link
Member Author

Choose a reason for hiding this comment

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

It’s not about a fixed number of Cells, rather it’s to stay a safe distance away from the convergence radius of the Taylor expansion.

Comment on lines +7 to +21
/** \brief calculate low integer powers base^exp
* \param[in] base base of power
*/
template<unsigned int exp> AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real pow (amrex::Real base) {
using namespace amrex::literals;
if constexpr (exp==0) {
return 1._rt;
} else if constexpr (exp==1) {
return base;
} else {
return pow<exp-1>(base) * base;
}
return 0._rt; //shut up compiler
}
Copy link
Member

Choose a reason for hiding this comment

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

std::pow is indeed much slower than multiplications for integer powers. @atmyers could that go to AMReX?

Comment on lines 149 to 153
* ``fields.extended_solve`` (`bool`) optional (default `0`)
Extends the area of the FFT Poisson solver to the ghost cells. This can reduce artefacts
originating from the boundary for long simulations.

* ``fields.open_boundary`` (`bool`) optional (default `0`)
Copy link
Member

Choose a reason for hiding this comment

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

Just for future extensibility, you could use the same syntax as in WarpX here: boundary.field_lo and boundary.field_hi
https://warpx.readthedocs.io/en/latest/usage/parameters.html#domain-boundary-conditions

Copy link
Member Author

Choose a reason for hiding this comment

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

No, because right now there is now functionality to have different boundary conditions for different sides of the box (and no use in hipace?). Moreover, as I understand it, the four combinations of fields.extended_solve = true/false and fields.open_boundary = true/false would require four differently name boundary conditions with that syntax.

@AlexanderSinn
Copy link
Member Author

CI test:
beam_in_vacuum

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.

Looks good to me, thanks, looking forward to testing it!

Copy link
Member

@SeverinDiederichs SeverinDiederichs left a comment

Choose a reason for hiding this comment

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

Great work! Thanks for this addition

@SeverinDiederichs SeverinDiederichs merged commit 93f01be into Hi-PACE:development Apr 5, 2022
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.

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants