Skip to content

Refactor Hartree_CG #387

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
jeanlucf22 merged 2 commits into from
Nov 26, 2025
Merged

Refactor Hartree_CG #387

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
884a21f
Template Hartree_CG on preconditioner data type
jeanlucf22 Nov 25, 2025
3440b41
Enable runtime precision selection for Hatree_CG preconditioner
jeanlucf22 Nov 26, 2025
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
22 changes: 13 additions & 9 deletions src/Control.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ Control::Control()
poisson_pc_nu1 = 2;
poisson_pc_nu2 = 2;
poisson_pc_nlev = 10;
poisson_pc_data_ = 32;
coloring_algo_ = 0;
maxDistanceAtomicInfo_ = 8.;
spread_factor = 2.;
Expand Down Expand Up @@ -331,7 +332,7 @@ void Control::sync(void)
if (onpe0 && verbose > 0)
(*MPIdata::sout) << "Control::sync()" << std::endl;
// pack
const short size_short_buffer = 91;
const short size_short_buffer = 92;
short* short_buffer = new short[size_short_buffer];
if (mype_ == 0)
{
Expand Down Expand Up @@ -421,6 +422,7 @@ void Control::sync(void)
short_buffer[88] = hartree_reset_;
short_buffer[89] = MD_last_step_;
short_buffer[90] = (short)static_cast<int>(poisson_lap_type_);
short_buffer[91] = poisson_pc_data_;
}
else
{
Expand Down Expand Up @@ -634,6 +636,7 @@ void Control::sync(void)
hartree_reset_ = short_buffer[88];
MD_last_step_ = short_buffer[89];
poisson_lap_type_ = static_cast<PoissonFDtype>(short_buffer[90]);
poisson_pc_data_ = short_buffer[91];

numst = int_buffer[0];
nel_ = int_buffer[1];
Expand Down Expand Up @@ -1406,14 +1409,15 @@ void Control::setOptions(const boost::program_options::variables_map& vm)
bool poisson_reset = vm["Poisson.reset"].as<bool>();
hartree_reset_ = poisson_reset ? 1 : 0;

poisson_pc_nu1 = vm["Poisson.nu1"].as<short>();
poisson_pc_nu2 = vm["Poisson.nu2"].as<short>();
vh_init = vm["Poisson.max_steps_initial"].as<short>();
vh_its = vm["Poisson.max_steps"].as<short>();
poisson_pc_nlev = vm["Poisson.max_levels"].as<short>();
rho0_ = vm["Poisson.rho0"].as<float>();
drho0_ = vm["Poisson.beta"].as<float>();
e0_ = vm["Poisson.e0"].as<float>();
poisson_pc_nu1 = vm["Poisson.nu1"].as<short>();
poisson_pc_nu2 = vm["Poisson.nu2"].as<short>();
vh_init = vm["Poisson.max_steps_initial"].as<short>();
vh_its = vm["Poisson.max_steps"].as<short>();
poisson_pc_nlev = vm["Poisson.max_levels"].as<short>();
rho0_ = vm["Poisson.rho0"].as<float>();
drho0_ = vm["Poisson.beta"].as<float>();
e0_ = vm["Poisson.e0"].as<float>();
poisson_pc_data_ = vm["Poisson.precond_precision"].as<short>();

str = vm["ProjectedMatrices.solver"].as<std::string>();
if (str.compare("short_sighted") == 0) short_sighted = 1;
Expand Down
5 changes: 5 additions & 0 deletions src/Control.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -403,6 +403,11 @@ class Control
short poisson_pc_nu2;
short poisson_pc_nlev;

/*!
* Poisson preconditioner precision (32 or 64)
*/
short poisson_pc_data_;

PoissonFDtype poisson_lap_type_;

short lap_type;
Expand Down
75 changes: 29 additions & 46 deletions src/Hartree_CG.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,13 +6,8 @@
// All rights reserved.
// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
// Please also read this link https://github.com/llnl/mgmol/LICENSE

#include <iomanip>
#include <iostream>
using namespace std;

#include "Control.h"
#include "Hartree_CG.h"
#include "Control.h"
#include "MultipoleExpansion.h"

#include "Laph2.h"
Expand All @@ -22,15 +17,16 @@ using namespace std;
#include "Laph6.h"
#include "Laph8.h"

// Timer Poisson::poisson_tm_("Poisson::poisson");
#include <iomanip>
#include <iostream>

template <class T>
void Hartree_CG<T>::solve(
const pb::GridFunc<RHODTYPE>& rho, const pb::GridFunc<RHODTYPE>& rhoc)
template <class OperType, typename ScalarType, typename PDataType>
void Hartree_CG<OperType, ScalarType, PDataType>::solve(
const pb::GridFunc<ScalarType>& rho, const pb::GridFunc<ScalarType>& rhoc)
{
PoissonInterface::poisson_tm_.start();

pb::GridFunc<RHODTYPE> work_rho(rho);
pb::GridFunc<ScalarType> work_rho(rho);
Control& ct = *(Control::instance());

// Keep in memory vh*rho before updating vh
Expand All @@ -48,7 +44,7 @@ void Hartree_CG<T>::solve(
if (Poisson::bc_[i] == 2) dim_mpol++;
//(*MPIdata::sout)<<"dim_mpol="<<dim_mpol<<endl;

pb::GridFunc<POTDTYPE> bc_func(
pb::GridFunc<ScalarType> bc_func(
Poisson::grid_, Poisson::bc_[0], Poisson::bc_[1], Poisson::bc_[2]);
if (dim_mpol > 0)
{
Expand All @@ -70,53 +66,40 @@ void Hartree_CG<T>::solve(
}
}

/* Check for uniform precision before calling poisson_solver.
* Downgrade or upgrade rhs (work_rho) to have precision of solution (vh_).
* Note that this could be done at the beginning of this function, but
* several operations involving rho might be done in lower precision
* (depending on POTDTYPE), which could affect accuracy. For now, we delay
* the switch until just before the solve call.
*/
// if(sizeof(POTDTYPE) != sizeof(RHODTYPE))
// {
/* solve with POTDTYPE precision */
pb::GridFunc<POTDTYPE> rhs(work_rho);
pb::GridFunc<ScalarType> rhs(work_rho);
rhs *= (4. * M_PI);
poisson_solver_->solve(*Poisson::vh_, rhs);
// }
// else
// {
// poisson_solver_->solve(*Poisson::vh_, work_rho);
// }

double residual_reduction = poisson_solver_->getResidualReduction();
double final_residual = poisson_solver_->getFinalResidual();

const double residual_reduction = poisson_solver_->getResidualReduction();
const double final_residual = poisson_solver_->getFinalResidual();
const bool large_residual
= (residual_reduction > 1.e-3 || final_residual > 1.e-3);

if (onpe0 && (large_residual || ct.verbose > 1))
(*MPIdata::sout) << setprecision(2) << scientific
(*MPIdata::sout) << std::setprecision(2) << std::scientific
<< "Hartree_CG: residual reduction = "
<< residual_reduction
<< ", final residual = " << final_residual << endl;
<< ", final residual = " << final_residual
<< std::endl;

Poisson::Int_vhrho_ = vel * Poisson::vh_->gdot(rho);
Poisson::Int_vhrhoc_ = vel * Poisson::vh_->gdot(rhoc);

PoissonInterface::poisson_tm_.stop();

assert(residual_reduction == residual_reduction);
assert(!std::isnan(residual_reduction));
}

template class Hartree_CG<pb::Laph2<POTDTYPE>>;
// template class Hartree_CG<pb::Laph2<float> >;
template class Hartree_CG<pb::Laph4<POTDTYPE>>;
// template class Hartree_CG<pb::Laph4<float> >;
template class Hartree_CG<pb::Laph6<POTDTYPE>>;
// template class Hartree_CG<pb::Laph6<float> >;
template class Hartree_CG<pb::Laph8<POTDTYPE>>;
// template class Hartree_CG<pb::Laph8<float> >;
template class Hartree_CG<pb::Laph4M<POTDTYPE>>;
// template class Hartree_CG<pb::Laph4M<float> >;
template class Hartree_CG<pb::Laph4MP<POTDTYPE>>;
// template class Hartree_CG<pb::Laph4MP<float> >;
template class Hartree_CG<pb::Laph2<double>, double, float>;
template class Hartree_CG<pb::Laph4<double>, double, float>;
template class Hartree_CG<pb::Laph6<double>, double, float>;
template class Hartree_CG<pb::Laph8<double>, double, float>;
template class Hartree_CG<pb::Laph4M<double>, double, float>;
template class Hartree_CG<pb::Laph4MP<double>, double, float>;

template class Hartree_CG<pb::Laph2<double>, double, double>;
template class Hartree_CG<pb::Laph4<double>, double, double>;
template class Hartree_CG<pb::Laph6<double>, double, double>;
template class Hartree_CG<pb::Laph8<double>, double, double>;
template class Hartree_CG<pb::Laph4M<double>, double, double>;
template class Hartree_CG<pb::Laph4MP<double>, double, double>;
23 changes: 14 additions & 9 deletions src/Hartree_CG.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,28 +7,33 @@
// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
// Please also read this link https://github.com/llnl/mgmol/LICENSE

#ifndef _HARTREE_CG_H_
#define _HARTREE_CG_H_
#ifndef MGMOL_HARTREE_CG_H
#define MGMOL_HARTREE_CG_H

#include "PCGSolver.h"
#include "Poisson.h"

template <class T>
#include <memory>

template <class OperatorType, typename ScalarType, typename PDataType>
class Hartree_CG : public Poisson
{
private:
PCGSolver<T, POTDTYPE>* poisson_solver_;
std::shared_ptr<PCGSolver<OperatorType, ScalarType, PDataType>>
poisson_solver_;

public:
// Constructor
Hartree_CG(const pb::Grid& grid, const short bc[3]) : Poisson(grid, bc)
{
T oper(Poisson::grid_);
poisson_solver_ = new PCGSolver<T, POTDTYPE>(oper, bc[0], bc[1], bc[2]);
OperatorType oper(Poisson::grid_);
poisson_solver_
= std::make_shared<PCGSolver<OperatorType, ScalarType, PDataType>>(
oper, bc[0], bc[1], bc[2]);
};

// Destructor
~Hartree_CG() override { delete poisson_solver_; }
~Hartree_CG() override {}

void setup(const short nu1, const short nu2, const short max_sweeps,
const double tol, const short max_nlevels,
Expand All @@ -38,8 +43,8 @@ class Hartree_CG : public Poisson
poisson_solver_->setup(nu1, nu2, max_sweeps, tol, max_nlevels);
}

void solve(const pb::GridFunc<RHODTYPE>& rho,
const pb::GridFunc<RHODTYPE>& rhoc) override;
void solve(const pb::GridFunc<ScalarType>& rho,
const pb::GridFunc<ScalarType>& rhoc) override;
};

#endif
Loading