One-dimensional Helmholtz kernel

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
   )
 

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 <https://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to MRCPP, see:
+ * <https://mrcpp.readthedocs.io/>
+ */
+
+#include <numeric>
+
+#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 <int D>
+Slater<D>::Slater(double a, double c, const Coord<D> &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 <int D> double Slater<D>::evalf(const Coord<D> &r) const {
+    auto dist = math_utils::calc_distance<D>(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

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 <https://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to MRCPP, see:
+ * <https://mrcpp.readthedocs.io/>
+ */
+
+/**
+ *
+ *  Base class for Slater type functions
+ */
+
+#pragma once
+
+#include <Eigen/Core>
+#include <cmath>
+#include <iostream>
+#include <memory>
+
+#include "MRCPP/mrcpp_declarations.h"
+#include "RepresentableFunction.h"
+
+namespace mrcpp {
+
+template <int D> class Slater : public RepresentableFunction<D> {
+public:
+    Slater(double a, double c, const Coord<D> &r);
+    Slater<D> &operator=(const Slater<D> &gp) = delete;
+    ~Slater() = default;
+
+    double evalf(const Coord<D> &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<double, D> &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<double, D> &r) { this->pos = r; }
+
+    friend std::ostream &operator<<(std::ostream &o, const Slater<D> &slater) { return slater.print(o); }
+
+protected:
+    double coef;                 /**< constant factor */
+    double alpha;                /**< exponent  */
+    Coord<D> 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

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 <https://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to MRCPP, see:
+ * <https://mrcpp.readthedocs.io/>
+ */
+
+#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

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 <https://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to MRCPP, see:
+ * <https://mrcpp.readthedocs.io/>
+ */
+
+#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
+

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)