[pull] master from MRChemSoft:master

api/MWOperators

@@ -38,3 +38,4 @@
 #include "operators/BSOperator.h"
 
 #include "treebuilders/apply.h"
+#include "treebuilders/DerivativeCalculator.h"

src/functions/RepresentableFunction.h

@@ -31,11 +31,13 @@
 
 #pragma once
 
+#include <Eigen/Core>
 #include <iostream>
 #include <vector>
 
 #include "MRCPP/constants.h"
 #include "MRCPP/mrcpp_declarations.h"
+#include "trees/NodeIndex.h"
 
 namespace mrcpp {
 
@@ -74,4 +76,15 @@ template <int D> class RepresentableFunction {
     virtual bool isZeroOnInterval(const double *a, const double *b) const { return false; }
 };
 
+/*
+ * Same as RepresentableFunction, but output a matrix of values
+ * for all points in a node, given its NodeIndex.
+ *
+ */
+class RepresentableFunction_M {
+public:
+    RepresentableFunction_M() {}
+    virtual Eigen::MatrixXd evalf(mrcpp::NodeIndex<3> nIdx) const = 0;
+};
+
 } // namespace mrcpp

src/treebuilders/DerivativeCalculator.cpp

@@ -86,6 +86,34 @@
     Printer::setPrecision(oldprec);
 }
 
+template <int D> void DerivativeCalculator<D>::calcNode(MWNode<D> &inpNode, MWNode<D> &outNode) {
+    //if (this->oper->getMaxBandWidth() > 1) MSG_ABORT("Only implemented for zero bw");
+    outNode.zeroCoefs();
+    int nComp = (1 << D);
+    double tmpCoefs[outNode.getNCoefs()];
+    OperatorState<D> os(outNode, tmpCoefs);
+
+    os.setFNode(inpNode);
+    os.setFIndex(inpNode.nodeIndex);
+    for (int ft = 0; ft < nComp; ft++) {
+        double fNorm = inpNode.getComponentNorm(ft);
+        if (fNorm < MachineZero) { continue; } // Could this be used outside the loop?
+        os.setFComponent(ft);
+        for (int gt = 0; gt < nComp; gt++) {
+            os.setGComponent(gt);
+            applyOperator_bw0(os);
+       }
+    }
+   // Multiply appropriate scaling factor. TODO: Could be included elsewhere
+    const double scaling_factor =
+        1.0/std::pow(outNode.getMWTree().getMRA().getWorldBox().getScalingFactor(this->applyDir), oper->getOrder());
+    if(abs(scaling_factor-1.0)>MachineZero){
+        for (int i = 0; i < outNode.getNCoefs(); i++) outNode.getCoefs()[i] *= scaling_factor;
+    }
+    outNode.calcNorms(); //TODO:required? norms are not used for now
+}
+
+
 template <int D> void DerivativeCalculator<D>::calcNode(MWNode<D> &gNode) {
     gNode.zeroCoefs();
 
@@ -101,6 +129,7 @@
     this->band_t[mrcpp_get_thread_num()].stop();
 
     this->calc_t[mrcpp_get_thread_num()].resume();
+
     for (int n = 0; n < fBand.size(); n++) {
         MWNode<D> &fNode = *fBand[n];
         NodeIndex<D> &fIdx = idx_band[n];
@@ -152,6 +181,39 @@
     return band;
 }
 
+/** Apply a single operator component (term) to a single f-node assuming zero bandwidth */
+template <int D> void DerivativeCalculator<D>::applyOperator_bw0(OperatorState<D> &os) {
+    //cout<<" applyOperator "<<endl;
+    MWNode<D> &gNode = *os.gNode;
+    MWNode<D> &fNode = *os.fNode;
+    const NodeIndex<D> &fIdx = *os.fIdx;
+    const NodeIndex<D> &gIdx = gNode.getNodeIndex();
+    int depth = gNode.getDepth();
+
+    double oNorm = 1.0;
+    double **oData = os.getOperData();
+
+    for (int d = 0; d < D; d++) {
+        const OperatorTree &oTree = this->oper->getComponent(0, d);
+        const OperatorNode &oNode = oTree.getNode(depth, 0);
+        int oIdx = os.getOperIndex(d);
+        if (this->applyDir == d) {
+            oData[d] = const_cast<double *>(oNode.getCoefs()) + oIdx * os.kp1_2;
+        } else {
+            if (oIdx == 0 or oIdx == 3) {
+                // This will activate the identity operator in direction i
+                oData[d] = nullptr;
+            } else {
+                // This means that we are in a zero part of the identity operator
+                return;
+            }
+        }
+    }
+    this->operStat.incrementFNodeCounters(fNode, os.ft, os.gt);
+    tensorApplyOperComp(os);
+}
+
+
 /** Apply a single operator component (term) to a single f-node. Whether the
 operator actualy is applied is determined by a screening threshold. */
 template <int D> void DerivativeCalculator<D>::applyOperator(OperatorState<D> &os) {
@@ -197,8 +259,8 @@
     tensorApplyOperComp(os);
 }
 
-/** Perorm the required linear algebra operations in order to apply an
-operator component to a f-node in a n-dimensional tesor space. */
+/** Perform the required linear algebra operations in order to apply an
+operator component to a f-node in a n-dimensional tensor space. */
 template <int D> void DerivativeCalculator<D>::tensorApplyOperComp(OperatorState<D> &os) {
     double **aux = os.getAuxData();
     double **oData = os.getOperData();

src/treebuilders/DerivativeCalculator.h

@@ -36,6 +36,7 @@ template <int D> class DerivativeCalculator final : public TreeCalculator<D> {
     ~DerivativeCalculator() override;
 
     MWNodeVector<D> *getInitialWorkVector(MWTree<D> &tree) const override;
+    void calcNode(MWNode<D> &fNode, MWNode<D> &gNode);
 
 private:
     int applyDir;
@@ -61,6 +62,7 @@ template <int D> class DerivativeCalculator final : public TreeCalculator<D> {
     }
 
     void applyOperator(OperatorState<D> &os);
+    void applyOperator_bw0(OperatorState<D> &os);
     void tensorApplyOperComp(OperatorState<D> &os);
 };
 

src/treebuilders/apply.cpp

@@ -256,6 +256,7 @@ template <int D> void apply_far_field(double prec, FunctionTree<D> &out, Convolu
 template <int D> void apply_near_field(double prec, FunctionTree<D> &out, ConvolutionOperator<D> &oper, FunctionTree<D> &inp, int maxIter, bool absPrec) {
     apply_on_unit_cell<D>(true, prec, out, oper, inp, maxIter, absPrec);
 }
+
 /** @brief Application of MW derivative operator
  *
  * @param[out] out: Output function to be built
@@ -274,7 +275,6 @@ template <int D> void apply_near_field(double prec, FunctionTree<D> &out, Convol
  */
 template <int D> void apply(FunctionTree<D> &out, DerivativeOperator<D> &oper, FunctionTree<D> &inp, int dir) {
     if (out.getMRA() != inp.getMRA()) MSG_ABORT("Incompatible MRA");
-
     TreeBuilder<D> builder;
     int maxScale = out.getMRA().getMaxScale();
 

src/trees/FunctionNode.cpp

@@ -130,7 +130,6 @@
     int qOrder = this->getKp1();
     getQuadratureCache(qc);
     const VectorXd &weights = qc.getWeights(qOrder);
-
     double sqWeights[qOrder];
     for (int i = 0; i < qOrder; i++) sqWeights[i] = std::sqrt(weights[i]);
 
@@ -158,6 +157,52 @@
     return two_n * sum;
 }
 
+/** Function integration, Interpolating basis.
+ *
+ * Integrates the function represented on the node on the full support of the
+ * node. A bit more involved than in the Legendre basis, as is requires some
+ * coupling of quadrature weights. */
+template <int D> double FunctionNode<D>::integrateValues() const {
+    int qOrder = this->getKp1();
+    getQuadratureCache(qc);
+    const VectorXd &weights = qc.getWeights(qOrder);
+    VectorXd coefs;
+    this->getCoefs(coefs);
+    int ncoefs = coefs.size();
+    int ncoefChild = ncoefs/(1<<D);
+    double cc[ncoefChild];
+    // factorize out the children
+    for (int i = 0; i < ncoefChild; i++)cc[i]=coefs[i];
+    for (int j = 1; j < (1<<D); j++) for (int i = 0; i < ncoefChild; i++)cc[i]+=coefs[j*ncoefChild+i];
+
+    int nc = 0;
+    double sum = 0.0;
+    if (D > 3) MSG_ABORT("Not Implemented")
+    else if (D == 3) {
+        for (int i = 0; i < qOrder; i++) {
+            double sumj = 0.0;
+            for (int j = 0; j < qOrder; j++) {
+                double sumk = 0.0;
+                for (int k = 0; k < qOrder; k++) sumk += cc[nc++] * weights[k];
+                sumj += sumk * weights[j];
+            }
+            sum += sumj * weights[i];
+        }
+    } else if (D==2) {
+        for (int j = 0; j < qOrder; j++) {
+                double sumk = 0.0;
+                for (int k = 0; k < qOrder; k++) sumk += cc[nc++] * weights[k];
+                sum += sumk * weights[j];
+        }
+    } else if (D==1) for (int k = 0; k < qOrder; k++) sum += cc[nc++] * weights[k];
+
+    int n = D * (this->getScale() + 1) ; // NB: one extra scale
+    int two_n = (1<<abs(n)); // 2**n;
+    if(n>0)sum/=two_n;
+    else sum*=two_n;
+    return sum;
+}
+
 template <int D> void FunctionNode<D>::setValues(const VectorXd &vec) {
     this->zeroCoefs();
     this->setCoefBlock(0, vec.size(), vec.data());
@@ -319,6 +364,9 @@
     auto &ftree = this->getFuncTree();
     if (this->isGenNode()) {
         ftree.getGenNodeAllocator().dealloc(sIdx);
+        // for GenNodes we clean unused chunks carefully, as they can become
+        // very large and occupy space long after used.
+        ftree.getGenNodeAllocator().deleteUnusedChunks();
     } else {
         ftree.decrementNodeCount(this->getScale());
         ftree.getNodeAllocator().dealloc(sIdx);

src/trees/FunctionNode.h

@@ -77,6 +77,7 @@ template <int D> class FunctionNode final : public MWNode<D> {
 
     double integrateLegendre() const;
     double integrateInterpolating() const;
+    double integrateValues() const;
 };
 
 template <int D> double dot_scaling(const FunctionNode<D> &bra, const FunctionNode<D> &ket);

src/trees/FunctionTree.cpp

@@ -54,7 +54,7 @@
 FunctionTree<D>::FunctionTree(const MultiResolutionAnalysis<D> &mra, SharedMemory *sh_mem, const std::string &name)
         : MWTree<D>(mra, name)
         , RepresentableFunction<D>(mra.getWorldBox().getLowerBounds().data(), mra.getWorldBox().getUpperBounds().data()) {
-    int nodesPerChunk = 64;
+    int nodesPerChunk = 2048; // Large chunks are required for not leading to memory fragmentation (32 MB on "Betzy" 2023)
     int coefsGenNodes = this->getKp1_d();
     int coefsRegNodes = this->getTDim() * this->getKp1_d();
     this->nodeAllocator_p = std::make_unique<NodeAllocator<D>>(this, sh_mem, coefsRegNodes, nodesPerChunk);
@@ -186,6 +186,42 @@
     return jacobian * result;
 }
 
+
+/** @returns Integral of a representable function over the grid given by the tree */
+template <> double FunctionTree<3>::integrateEndNodes(RepresentableFunction_M &f) {
+    //traverse tree, and treat end nodes only
+    std::vector<FunctionNode<3> *> stack;   // node from this
+    for (int i = 0; i < this->getRootBox().size(); i++) stack.push_back(&(this->getRootFuncNode(i)));
+    int basis = getMRA().getScalingBasis().getScalingType();
+    double result = 0.0;
+    while (stack.size() > 0) {
+        FunctionNode<3> *Node = stack.back();
+        stack.pop_back();
+        if (Node->getNChildren() > 0) {
+            for (int i = 0; i < Node->getNChildren(); i++) stack.push_back(&(Node->getFuncChild(i)));
+        } else {
+            //end nodes
+            Eigen::MatrixXd fmat = f.evalf(Node->nodeIndex);
+            double *coefs = Node->getCoefs(); // save position of coeff, but do not use them!
+            // The data in fmat is not organized so that two consecutive points are stored after each other in memory, so needs to copy before mwtransform, cannot use memory adress directly.
+            int nc=fmat.cols();
+            double cc[nc];
+            for (int i = 0; i < nc; i++)cc[i]=fmat(0,i);
+            Node->attachCoefs(cc);
+            result += Node->integrateValues();
+            Node->attachCoefs(coefs); // put back original coeff
+        }
+    }
+
+    // Handle potential scaling
+    auto scaling_factor = this->getMRA().getWorldBox().getScalingFactors();
+    auto jacobian = 1.0;
+    for (const auto &sf_i : scaling_factor) jacobian *= std::sqrt(sf_i);
+    // Square root of scaling factor in each diection. The seemingly missing
+    // multiplication by the square root of sf_i is included in the basis
+    return jacobian * result;
+}
+
 /** @returns Function value in a point, out of bounds returns zero
  *
  * @param[in] r: Cartesian coordinate

src/trees/FunctionTree.h

@@ -62,6 +62,7 @@ template <int D> class FunctionTree final : public MWTree<D>, public Representab
     ~FunctionTree() override;
 
     double integrate() const;
+    double integrateEndNodes(RepresentableFunction_M &f);
     double evalf_precise(const Coord<D> &r);
     double evalf(const Coord<D> &r) const override;
 

src/trees/MWNode.cpp

@@ -103,7 +103,7 @@ MWNode<D>::MWNode(MWTree<D> *tree, int rIdx)
 /** @brief MWNode constructor.
  *
  * @param[in] parent: parent node
- * @param[in] cIdx: child index of the current node 
+ * @param[in] cIdx: child index of the current node
  *
  * @details Constructor for leaf nodes. It requires the corresponding
  * parent and an integer to identify the correct child.
@@ -131,7 +131,7 @@ MWNode<D>::MWNode(MWNode<D> *parent, int cIdx)
  * the tree.
  */
 template <int D>
-MWNode<D>::MWNode(const MWNode<D> &node, bool allocCoef)
+MWNode<D>::MWNode(const MWNode<D> &node, bool allocCoef, bool SetCoef)
         : tree(node.tree)
         , parent(nullptr)
         , nodeIndex(node.nodeIndex)
@@ -141,7 +141,7 @@ MWNode<D>::MWNode(const MWNode<D> &node, bool allocCoef)
     setIsLooseNode();
     if (allocCoef) {
         allocCoefs(this->getTDim(), this->getKp1_d());
-        if (node.hasCoefs()) {
+        if (node.hasCoefs() and SetCoef) {
             setCoefBlock(0, node.getNCoefs(), node.getCoefs());
             if (this->getNCoefs() > node.getNCoefs()) {
                 for (int i = node.getNCoefs(); i < this->getNCoefs(); i++) this->coefs[i] = 0.0;
@@ -236,7 +236,7 @@ template <int D> void MWNode<D>::getCoefs(Eigen::VectorXd &c) const {
     c = VectorXd::Map(this->coefs, this->n_coefs);
 }
 
-/** @brief sets all MW coefficients and the norms to zero 
+/** @brief sets all MW coefficients and the norms to zero
  *
  */
 template <int D> void MWNode<D>::zeroCoefs() {
@@ -388,7 +388,7 @@ template <int D> void MWNode<D>::giveParentCoefs(bool overwrite) {
 }
 
 /** @brief Copy scaling coefficients from children to parent
- * 
+ *
  * @details Takes the scaling coefficients of the children and stores
  * them consecutively in the corresponding block of the parent,
  * following the usual bitwise notation.
@@ -404,7 +404,7 @@ template <int D> void MWNode<D>::copyCoefsFromChildren() {
 }
 
 /** @brief Generates scaling cofficients of children
- * 
+ *
  * @details If the node is a leafNode, it takes the scaling&wavelet
  * coefficients of the parent and it generates the scaling
  * coefficients for the children
@@ -835,7 +835,7 @@ template <int D> void MWNode<D>::getPrimitiveQuadPts(MatrixXd &pts) const {
  * the set of quadrature points becomes \f$ x^\alpha_i = 2^{-n-1} (x_i + 2 l^\alpha + t^\alpha) \f$, where \f$ t^\alpha =
  * 0,1 \f$. By taking all possible \f$(k+1)^d\combinations \f$, they will
  * then define a d-dimensional grid of quadrature points for the child
- * nodes. 
+ * nodes.
  *
  */
 template <int D> void MWNode<D>::getPrimitiveChildPts(MatrixXd &pts) const {
@@ -879,7 +879,7 @@ template <int D> void MWNode<D>::getExpandedQuadPts(Eigen::MatrixXd &pts) const
 
 /** @brief Returns the quadrature points in a given node
  *
- * @param[in,out] pts: expanded quadrature points in a \f$ d \times 
+ * @param[in,out] pts: expanded quadrature points in a \f$ d \times
  * 2^d(k+1)^d \f$ matrix form.
  *
  * @details The primitive quadrature points of the children are used to obtain a