diff --git a/src/functions/CMakeLists.txt b/src/functions/CMakeLists.txt index c55092d81..2e4cd879d 100644 --- a/src/functions/CMakeLists.txt +++ b/src/functions/CMakeLists.txt @@ -10,6 +10,7 @@ target_sources(mrcpp ${CMAKE_CURRENT_SOURCE_DIR}/GaussFunc.cpp ${CMAKE_CURRENT_SOURCE_DIR}/GaussPoly.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Gaussian.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/Slater.cpp ${CMAKE_CURRENT_SOURCE_DIR}/special_functions.cpp ) @@ -27,6 +28,7 @@ list(APPEND ${_dirname}_h ${CMAKE_CURRENT_SOURCE_DIR}/GaussFunc.h ${CMAKE_CURRENT_SOURCE_DIR}/GaussPoly.h ${CMAKE_CURRENT_SOURCE_DIR}/Gaussian.h + ${CMAKE_CURRENT_SOURCE_DIR}/Slater.h ${CMAKE_CURRENT_SOURCE_DIR}/special_functions.h ) diff --git a/src/functions/Slater.cpp b/src/functions/Slater.cpp new file mode 100644 index 000000000..e78e22bba --- /dev/null +++ b/src/functions/Slater.cpp @@ -0,0 +1,74 @@ +/* + * MRCPP, a numerical library based on multiresolution analysis and + * the multiwavelet basis which provide low-scaling algorithms as well as + * rigorous error control in numerical computations. + * Copyright (C) 2021 Stig Rune Jensen, Jonas Juselius, Luca Frediani and contributors. + * + * This file is part of MRCPP. + * + * MRCPP is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRCPP is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRCPP. If not, see . + * + * For information on the complete list of contributors to MRCPP, see: + * + */ + +#include + +#include "Slater.h" +#include "function_utils.h" +#include "trees/NodeIndex.h" +#include "utils/Printer.h" +#include "utils/details.h" +#include "utils/math_utils.h" + +using namespace Eigen; + +namespace mrcpp { + +template +Slater::Slater(double a, double c, const Coord &r) + : coef(c) + , alpha(a) + , pos(r) { +} + +template<> double Slater<1>::calcSquareNorm() const { + double c2 = this->coef * this->coef; + double square_norm = c2 / this->alpha; + return square_norm; +} + +template<> double Slater<2>::calcSquareNorm() const { + double c2 = this->coef * this->coef; + double k2 = this->alpha * this->alpha; + return 0.5 * pi * c2 / k2; +} + +template<> double Slater<3>::calcSquareNorm() const { + double c2 = this->coef * this->coef; + double k3 = this->alpha * this->alpha * this->alpha; + return pi * c2 / k3; +} + +template double Slater::evalf(const Coord &r) const { + auto dist = math_utils::calc_distance(r , pos); + double exp_val = std::exp(-this->alpha * dist); + return this->coef * exp_val; +} + +template class Slater<1>; +template class Slater<2>; +template class Slater<3>; + +} // namespace mrcpp diff --git a/src/functions/Slater.h b/src/functions/Slater.h new file mode 100644 index 000000000..081833453 --- /dev/null +++ b/src/functions/Slater.h @@ -0,0 +1,86 @@ +/* + * MRCPP, a numerical library based on multiresolution analysis and + * the multiwavelet basis which provide low-scaling algorithms as well as + * rigorous error control in numerical computations. + * Copyright (C) 2021 Stig Rune Jensen, Jonas Juselius, Luca Frediani and contributors. + * + * This file is part of MRCPP. + * + * MRCPP is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRCPP is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRCPP. If not, see . + * + * For information on the complete list of contributors to MRCPP, see: + * + */ + +/** + * + * Base class for Slater type functions + */ + +#pragma once + +#include +#include +#include +#include + +#include "MRCPP/mrcpp_declarations.h" +#include "RepresentableFunction.h" + +namespace mrcpp { + +template class Slater : public RepresentableFunction { +public: + Slater(double a, double c, const Coord &r); + Slater &operator=(const Slater &gp) = delete; + ~Slater() = default; + + double evalf(const Coord &r) const; + + double calcSquareNorm() const; + + /** @brief Rescale function by its norm \f$ ||f||^{-1} \f$ */ + void normalize() { + double norm = std::sqrt(calcSquareNorm()); + multConstInPlace(1.0 / norm); + } + void multConstInPlace(double c) { this->coef *= c; } + void operator*=(double c) { multConstInPlace(c); } + + // bool getScreen() const { return screen; } + // bool checkScreen(int n, const int *l) const; + + double getCoef() const { return coef; } + double getAlpha() const { return alpha; } + const std::array &getPos() const { return pos; } + + // void setScreen(bool _screen) { this->screen = _screen; } + void setCoef(double cf) { this->coef = cf; } + void setAlpha(double _alpha) { this->alpha = _alpha; } + void setPos(const std::array &r) { this->pos = r; } + + friend std::ostream &operator<<(std::ostream &o, const Slater &slater) { return slater.print(o); } + +protected: + double coef; /**< constant factor */ + double alpha; /**< exponent */ + Coord pos; /**< center */ + + bool isVisibleAtScale(int scale, int nQuadPts) const {return true;} + bool isZeroOnInterval(const double *a, const double *b) const {return false;} + + // virtual std::ostream &print(std::ostream &o) const = 0; +}; + +} // namespace mrcpp diff --git a/src/operators/HelmholtzOperator1d.cpp b/src/operators/HelmholtzOperator1d.cpp new file mode 100644 index 000000000..3225cefe8 --- /dev/null +++ b/src/operators/HelmholtzOperator1d.cpp @@ -0,0 +1,80 @@ +/* + * MRCPP, a numerical library based on multiresolution analysis and + * the multiwavelet basis which provide low-scaling algorithms as well as + * rigorous error control in numerical computations. + * Copyright (C) 2021 Stig Rune Jensen, Jonas Juselius, Luca Frediani and contributors. + * + * This file is part of MRCPP. + * + * MRCPP is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRCPP is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRCPP. If not, see . + * + * For information on the complete list of contributors to MRCPP, see: + * + */ + +#include "HelmholtzOperator1d.h" +#include "HelmholtzKernel.h" +#include "utils/Printer.h" + +namespace mrcpp { + +/** @returns New HelmholtzOperator object + * @param[in] mra: Which MRA the operator is defined + * @param[in] m: Exponential parameter of the operator + * @param[in] pr: Build precision, closeness to exp(-mu*r)/r + * @details This will construct a gaussian expansion to approximate + * exp(-mu*r)/r, and project each term into a one-dimensional MW operator. + * Subsequent application of this operator will apply each of the terms to + * the input function in all Cartesian directions. + */ +HelmholtzOperator::HelmholtzOperator(const MultiResolutionAnalysis<3> &mra, double mu, double prec) + : ConvolutionOperator<3>(mra) { + int oldlevel = Printer::setPrintLevel(0); + + this->setBuildPrec(prec); + double o_prec = prec; + double k_prec = prec / 10.0; + double r_min = this->MRA.calcMinDistance(k_prec); + double r_max = this->MRA.calcMaxDistance(); + + HelmholtzKernel kernel(mu, k_prec, r_min, r_max); + initialize(kernel, k_prec, o_prec); + this->initOperExp(kernel.size()); + + Printer::setPrintLevel(oldlevel); +} + +HelmholtzOperator::HelmholtzOperator(const MultiResolutionAnalysis<3> &mra, double mu, double prec, int root, int reach) + : ConvolutionOperator<3>(mra, root, reach) { + int oldlevel = Printer::setPrintLevel(0); + + this->setBuildPrec(prec); + double o_prec = prec; + double k_prec = prec / 100.0; + double r_min = this->MRA.calcMinDistance(k_prec); + double r_max = this->MRA.calcMaxDistance(); + + // Adjust r_max for periodic world + auto rel_root = this->oper_root - this->MRA.getRootScale(); + r_max *= std::pow(2.0, -rel_root); + r_max *= (2.0 * this->oper_reach) + 1.0; + + HelmholtzKernel kernel(mu, k_prec, r_min, r_max); + initialize(kernel, k_prec, o_prec); + this->initOperExp(kernel.size()); + + Printer::setPrintLevel(oldlevel); +} + +} // namespace mrcpp diff --git a/src/operators/HelmholtzOperator1d.h b/src/operators/HelmholtzOperator1d.h new file mode 100644 index 000000000..c052e9326 --- /dev/null +++ b/src/operators/HelmholtzOperator1d.h @@ -0,0 +1,51 @@ +/* + * MRCPP, a numerical library based on multiresolution analysis and + * the multiwavelet basis which provide low-scaling algorithms as well as + * rigorous error control in numerical computations. + * Copyright (C) 2021 Stig Rune Jensen, Jonas Juselius, Luca Frediani and contributors. + * + * This file is part of MRCPP. + * + * MRCPP is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRCPP is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRCPP. If not, see . + * + * For information on the complete list of contributors to MRCPP, see: + * + */ + +#pragma once + +#include "ConvolutionOperator.h" + +namespace mrcpp { + +/** @class HelmholtzOperator1d + * + * @brief Convolution with the Helmholtz Green's function kernel for the one-dimensional case + * + * @details The 1d Helmholtz kernel is a single Slater-type function + * + * \f$ H(r-r') = \frac{1}{2k} e^{-k}{|x-x'|} \f$ + */ + + +class HelmholtzOperator1d final : public ConvolutionOperator<3> { +public: + HelmholtzOperator1d(const MultiResolutionAnalysis<3> &mra, double k, double prec); + HelmholtzOperator1d(const MultiResolutionAnalysis<3> &mra, double k, double prec, int root, int reach = 1); + HelmholtzOperator1d(const HelmholtzOperator1d &oper) = delete; + HelmholtzOperator1d &operator=(const HelmholtzOperator1d &oper) = delete; +}; + +} // namespace mrcpp + diff --git a/tests/functions/CMakeLists.txt b/tests/functions/CMakeLists.txt index e322a7fd2..b6d10372d 100644 --- a/tests/functions/CMakeLists.txt +++ b/tests/functions/CMakeLists.txt @@ -2,12 +2,14 @@ target_sources(mrcpp-tests PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/legendre_poly.cpp ${CMAKE_CURRENT_SOURCE_DIR}/polynomial.cpp ${CMAKE_CURRENT_SOURCE_DIR}/gaussians.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/slater.cpp ${CMAKE_CURRENT_SOURCE_DIR}/periodify_gaussians.cpp ) -add_Catch_test(NAME legendre_poly LABELS legendre_poly) -add_Catch_test(NAME polynomials LABELS polynomials) -add_Catch_test(NAME gaussians LABELS gaussians) -add_Catch_test(NAME periodic_narrow_gaussian LABELS periodic_narrow_gaussian) -add_Catch_test(NAME periodic_wide_gaussian LABELS periodic_wide_gaussian) -add_Catch_test(NAME periodic_gaussExp LABELS periodic_gausExp) +add_Catch_test(NAME legendre_poly LABELS legendre_poly) +add_Catch_test(NAME polynomials LABELS polynomials) +add_Catch_test(NAME gaussians LABELS gaussians) +add_Catch_test(NAME slater LABELS slater) +add_Catch_test(NAME periodic_narrow_gaussian LABELS periodic_narrow_gaussian) +add_Catch_test(NAME periodic_wide_gaussian LABELS periodic_wide_gaussian) +add_Catch_test(NAME periodic_gaussExp LABELS periodic_gausExp) diff --git a/tests/functions/slater.cpp b/tests/functions/slater.cpp new file mode 100644 index 000000000..aa578b697 --- /dev/null +++ b/tests/functions/slater.cpp @@ -0,0 +1,126 @@ +/* + * MRCPP, a numerical library based on multiresolution analysis and + * the multiwavelet basis which provide low-scaling algorithms as well as + * rigorous error control in numerical computations. + * Copyright (C) 2021 Stig Rune Jensen, Jonas Juselius, Luca Frediani and contributors. + * + * This file is part of MRCPP. + * + * MRCPP is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRCPP is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRCPP. If not, see . + * + * For information on the complete list of contributors to MRCPP, see: + * + */ + +#include "catch2/catch_all.hpp" + +#include "factory_functions.h" + +#include "functions/Slater.h" +#include "treebuilders/grid.h" +#include "treebuilders/project.h" +#include "utils/math_utils.h" + +using namespace mrcpp; + +namespace gaussians { + +template void testSlaterFunction(); + +SCENARIO("Slater", "[slater]") { + GIVEN("Slater Function in 1D") { testSlaterFunction<1>(); } + GIVEN("Slater Function in 2D") { testSlaterFunction<2>(); } + GIVEN("Slater Function in 3D") { testSlaterFunction<3>(); } +} + +template void testSlaterFunction() { + const auto prec = 1.0e-4; + // Making normalized gaussian centered at the origin + const auto alpha = 2.; + const auto coeff = 3.; + double pos_data[3] = {0.1, 0.2, 0.3}; + double ref_data[3] = {-0.2, 0.5, 1.0}; + auto pos = mrcpp::details::convert_to_std_array(pos_data); + auto ref = mrcpp::details::convert_to_std_array(ref_data); + + auto slater = Slater(alpha, coeff, pos); + + // Making ref value + auto dist = math_utils::calc_distance(slater.getPos(), ref); + const auto ref_val = slater.getCoef() * std::exp(-slater.getAlpha() * std::abs(dist)); + + // Making MRA + const auto min_scale = -4; + auto corner = std::array{}; + auto boxes = std::array{}; + corner.fill(-1); + boxes.fill(2); + auto world = BoundingBox(min_scale, corner, boxes); + const auto basis = InterpolatingBasis(7); + auto MRA = MultiResolutionAnalysis(world, basis, 25); + + // Making a function tree and projects the slater fcn onto it + FunctionTree f_tree(MRA); + build_grid(f_tree, slater); + project(prec, f_tree, slater); + f_tree.calcSquareNorm(); + + WHEN("Slater function is projected") { + THEN("The Slater function can be evaluated") { REQUIRE(slater.evalf(ref) == Catch::Approx(ref_val)); } + THEN("The projected Slater function can be evaluated") { REQUIRE(f_tree.evalf(ref) == Catch::Approx(ref_val)); } + THEN("the square norm matches the analytical value") { REQUIRE(f_tree.getSquareNorm() == Catch::Approx(slater.calcSquareNorm())); } + } +} + + /* +SCENARIO("Slater_2D", "[slater_2d]") { + const auto D = 2; + const auto prec = 1.0e-4; + // Making normalized gaussian centered at the origin + const auto pos = Coord{0.0, 0.0}; + const auto alpha = 2.2; + const auto coeff = 3.3; + auto slater = Slater(alpha, coeff, pos); + + // Making ref value + const auto r_ref = std::array{.1, 0.2}; + double dist = math_utils::calc_distance(r_ref, slater.getPos()); + const auto ref_val = slater.getCoef() * std::exp(-slater.getAlpha() * std::abs(dist)); + + // Making MRA + const auto min_scale = -4; + const auto corner = std::array{-1,-1}; + const auto boxes = std::array{2,2}; + auto world = BoundingBox(min_scale, corner, boxes); + const auto basis = InterpolatingBasis(7); + auto MRA = MultiResolutionAnalysis(world, basis, 25); + + // Making a function tree and projects the slater fcn onto it + FunctionTree f_tree(MRA); + build_grid(f_tree, slater); + project(prec, f_tree, slater); + f_tree.calcSquareNorm(); + + WHEN("Slater function is projected") { + THEN("The Slater function can be evaluated") { REQUIRE(slater.evalf(r_ref) == Catch::Approx(ref_val)); } + + THEN("The projected Slater function can be evaluated") { REQUIRE(f_tree.evalf(r_ref) == Catch::Approx(ref_val)); } + + THEN("the square norm matches") { REQUIRE(f_tree.getSquareNorm() == Catch::Approx(slater.calcSquareNorm())); } + + } +} + */ + +} // namespace gaussians