Skip to content

Commit

Permalink
Merge branch 'free_function_poisson_solver' into 'master'
Browse files Browse the repository at this point in the history
Make the poisson solve a free function instead of using an empty class

See merge request npneq/inq!1149
  • Loading branch information
xavierandrade committed Oct 15, 2024
2 parents 1315504 + e8a27e1 commit 684a808
Show file tree
Hide file tree
Showing 5 changed files with 28 additions and 42 deletions.
6 changes: 2 additions & 4 deletions python/_pinq.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -136,10 +136,8 @@ struct calculator {
electrons_.emplace(systems::electrons(ions, els_));
ground_state::initial_guess(ions, *electrons_);
}

solvers::poisson poisson_solver;

auto ionic_long_range = poisson_solver(electrons_->atomic_pot().ionic_density(electrons_->kpin_states_comm(), electrons_->density_basis(), ions));

auto ionic_long_range = solvers::poisson::solve(electrons_->atomic_pot().ionic_density(electrons_->kpin_states_comm(), electrons_->density_basis(), ions));
auto ionic_short_range = electrons_->atomic_pot().local_potential(electrons_->kpin_states_comm(), electrons_->density_basis(), ions);
auto vion = operations::add(ionic_long_range, ionic_short_range);

Expand Down
3 changes: 1 addition & 2 deletions src/hamiltonian/exchange_operator.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,6 @@ namespace hamiltonian {
gpu::array<int, 1> kpoint_indices_;
std::optional<basis::field_set<basis::real_space, complex, parallel::arbitrary_partition>> orbitals_;
std::vector<states::orbital_set<basis::real_space, complex>> ace_orbitals_;
solvers::poisson poisson_solver_;
double exchange_coefficient_;
bool use_ace_;
singularity_correction sing_;
Expand Down Expand Up @@ -142,7 +141,7 @@ namespace hamiltonian {
});
}

poisson_solver_.in_place(rhoij, -phi.kpoint() + kpt[jj], sing_(idx[jj]));
solvers::poisson::in_place(rhoij, -phi.kpoint() + kpt[jj], sing_(idx[jj]));

{ CALI_CXX_MARK_SCOPE("exchange_operator::mulitplication");
gpu::run(nst, exxphi.basis().local_size(),
Expand Down
8 changes: 2 additions & 6 deletions src/hamiltonian/self_consistency.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -74,9 +74,7 @@ class self_consistency {

CALI_CXX_MARK_FUNCTION;

solvers::poisson poisson_solver;

auto ionic_long_range = poisson_solver(atomic_pot.ionic_density(comm, density_basis_, ions));
auto ionic_long_range = solvers::poisson::solve(atomic_pot.ionic_density(comm, density_basis_, ions));
auto ionic_short_range = atomic_pot.local_potential(comm, density_basis_, ions);
vion_ = operations::add(ionic_long_range, ionic_short_range);

Expand All @@ -97,8 +95,6 @@ class self_consistency {

energy.external(operations::integral_product(total_density, vion_));

solvers::poisson poisson_solver;

//IONIC POTENTIAL
auto vscalar = vion_;

Expand All @@ -117,7 +113,7 @@ class self_consistency {

// Hartree
if(theory_.hartree_potential()){
auto vhartree = poisson_solver(total_density);
auto vhartree = solvers::poisson::solve(total_density);
energy.hartree(0.5*operations::integral_product(total_density, vhartree));
operations::increment(vscalar, vhartree);
} else {
Expand Down
4 changes: 1 addition & 3 deletions src/observables/forces_stress.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -81,11 +81,9 @@ struct forces_stress {

{ CALI_CXX_MARK_SCOPE("forces_local");

solvers::poisson poisson_solver;

//the force from the local potential
for(int iatom = 0; iatom < ions.size(); iatom++){
auto ionic_long_range = poisson_solver(electrons.atomic_pot().ionic_density(electrons.states_comm(), electrons.density_basis(), ions, iatom));
auto ionic_long_range = solvers::poisson::solve(electrons.atomic_pot().ionic_density(electrons.states_comm(), electrons.density_basis(), ions, iatom));
auto ionic_short_range = electrons.atomic_pot().local_potential(electrons.states_comm(), electrons.density_basis(), ions, iatom);

auto force_cov = -gpu::run(gpu::reduce(electrons.density_basis().local_size()),
Expand Down
49 changes: 22 additions & 27 deletions src/solvers/poisson.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,8 @@ class poisson {

public:

poisson() = delete;

struct poisson_kernel_3d {
GPU_FUNCTION auto operator()(vector3<double, cartesian> gg, double const zeroterm) const {
auto g2 = norm(gg);
Expand Down Expand Up @@ -67,7 +69,7 @@ class poisson {
///////////////////////////////////////////////////////////////////////////////////////////////////

template <typename KernelType, typename FieldSetType>
void poisson_apply_kernel(KernelType const kernel, FieldSetType & density, vector3<double> const & gshift = {0.0, 0.0, 0.0}, double const zeroterm = 0.0) const {
static void poisson_apply_kernel(KernelType const kernel, FieldSetType & density, vector3<double> const & gshift = {0.0, 0.0, 0.0}, double const zeroterm = 0.0) {

static_assert(std::is_same_v<typename FieldSetType::basis_type, basis::fourier_space>, "Only makes sense in fourier_space");

Expand All @@ -87,7 +89,7 @@ class poisson {

private:

auto poisson_solve_3d(basis::field<basis::real_space, complex> const & density) const {
static auto poisson_solve_3d(basis::field<basis::real_space, complex> const & density) {

CALI_CXX_MARK_FUNCTION;

Expand All @@ -98,7 +100,7 @@ class poisson {

///////////////////////////////////////////////////////////////////////////////////////////////////

void poisson_solve_in_place_3d(basis::field_set<basis::real_space, complex> & density, vector3<double> const & gshift, double const zeroterm) const {
static void poisson_solve_in_place_3d(basis::field_set<basis::real_space, complex> & density, vector3<double> const & gshift, double const zeroterm) {

CALI_CXX_MARK_FUNCTION;

Expand All @@ -109,7 +111,7 @@ class poisson {

///////////////////////////////////////////////////////////////////////////////////////////////////

basis::field<basis::real_space, complex> poisson_solve_2d(basis::field<basis::real_space, complex> const & density) const {
static basis::field<basis::real_space, complex> poisson_solve_2d(basis::field<basis::real_space, complex> const & density) {

CALI_CXX_MARK_FUNCTION;

Expand All @@ -127,7 +129,7 @@ class poisson {

///////////////////////////////////////////////////////////////////////////////////////////////////

void poisson_solve_in_place_2d(basis::field_set<basis::real_space, complex> & density, vector3<double> const & gshift, double const zeroterm) const {
static void poisson_solve_in_place_2d(basis::field_set<basis::real_space, complex> & density, vector3<double> const & gshift, double const zeroterm) {

CALI_CXX_MARK_FUNCTION;

Expand All @@ -143,7 +145,7 @@ class poisson {

///////////////////////////////////////////////////////////////////////////////////////////////////

basis::field<basis::real_space, complex> poisson_solve_0d(basis::field<basis::real_space, complex> const & density) const {
static basis::field<basis::real_space, complex> poisson_solve_0d(basis::field<basis::real_space, complex> const & density) {

CALI_CXX_MARK_FUNCTION;

Expand All @@ -161,7 +163,7 @@ class poisson {

///////////////////////////////////////////////////////////////////////////////////////////////////

void poisson_solve_in_place_0d(basis::field_set<basis::real_space, complex> & density, vector3<double> const & gshift, double const zeroterm) const {
static void poisson_solve_in_place_0d(basis::field_set<basis::real_space, complex> & density, vector3<double> const & gshift, double const zeroterm) {

CALI_CXX_MARK_FUNCTION;

Expand All @@ -179,7 +181,7 @@ class poisson {

public:

auto operator()(const basis::field<basis::real_space, complex> & density) const {
static auto solve(const basis::field<basis::real_space, complex> & density) {

CALI_CXX_MARK_SCOPE("poisson(complex)");

Expand All @@ -195,7 +197,7 @@ class poisson {
///////////////////////////////////////////////////////////////////////////////////////////////////

template <typename Space = cartesian>
void in_place(basis::field_set<basis::real_space, complex> & density, vector3<double, Space> const & gshift = {0.0, 0.0, 0.0}, double const zeroterm = 0.0) const {
static void in_place(basis::field_set<basis::real_space, complex> & density, vector3<double, Space> const & gshift = {0.0, 0.0, 0.0}, double const zeroterm = 0.0) {

CALI_CXX_MARK_SCOPE("poisson(complex)");

Expand All @@ -212,11 +214,11 @@ class poisson {

///////////////////////////////////////////////////////////////////////////////////////////////////

basis::field<basis::real_space, double> operator()(const basis::field<basis::real_space, double> & density) const {
static basis::field<basis::real_space, double> solve(const basis::field<basis::real_space, double> & density) {

CALI_CXX_MARK_SCOPE("poisson(real)");

auto complex_potential = operator()(complex_field(density));
auto complex_potential = solve(complex_field(density));
return real_field(complex_potential);
}

Expand All @@ -234,8 +236,6 @@ class poisson {
#include <operations/integral.hpp>

TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {


using namespace inq;
using namespace inq::magnitude;
using namespace Catch::literals;
Expand All @@ -262,7 +262,6 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {

field<real_space, complex> density(rs);
field_set<real_space, complex> density_set(rs, nst);
solvers::poisson psolver;

SECTION("Point charge"){

Expand All @@ -284,8 +283,8 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {
}
}

auto potential = psolver(density);
psolver.in_place(density_set);
auto potential = solvers::poisson::solve(density);
solvers::poisson::in_place(density_set);

double sum[2] = {0.0, 0.0};
for(int ix = 0; ix < rs.local_sizes()[0]; ix++){
Expand Down Expand Up @@ -332,8 +331,8 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {
}
}

auto potential = psolver(density);
psolver.in_place(density_set);
auto potential = solvers::poisson::solve(density);
solvers::poisson::in_place(density_set);

double diff = 0.0;
for(int ix = 0; ix < rs.local_sizes()[0]; ix++){
Expand Down Expand Up @@ -372,7 +371,7 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {
}
}

auto rpotential = psolver(rdensity);
auto rpotential = solvers::poisson::solve(rdensity);

double diff = 0.0;
for(int ix = 0; ix < rs.local_sizes()[0]; ix++){
Expand All @@ -396,8 +395,6 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {
{
basis::real_space rs(systems::cell::cubic(8.0_b).finite(), /*spacing =*/ 0.09, comm);

solvers::poisson psolver;

SECTION("Grid finite"){

CHECK(rs.cell().periodicity() == 0);
Expand Down Expand Up @@ -430,8 +427,8 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {

CHECK(real(operations::integral(density)) == -1.0_a);

auto potential = psolver(density);
psolver.in_place(density_set);
auto potential = solvers::poisson::solve(density);
solvers::poisson::in_place(density_set);

for(int ix = 0; ix < rs.local_sizes()[0]; ix++){
for(int iy = 0; iy < rs.local_sizes()[1]; iy++){
Expand Down Expand Up @@ -476,8 +473,6 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {

basis::real_space rs(systems::cell::orthorhombic(6.0_b, 6.0_b, 9.0_b).periodicity(2), /*spacing =*/ 0.12, comm);

solvers::poisson psolver;

CHECK(rs.cell().periodicity() == 2);

CHECK(rs.sizes()[0] == 50);
Expand All @@ -504,8 +499,8 @@ TEST_CASE(INQ_TEST_FILE, INQ_TEST_TAG) {

CHECK(real(operations::integral(density)) == -1.0_a);

auto potential = psolver(density);
psolver.in_place(density_set);
auto potential = solvers::poisson::solve(density);
solvers::poisson::in_place(density_set);

auto & part = potential.basis().part();

Expand Down

0 comments on commit 684a808

Please sign in to comment.