Merge branch 'master' into one-corner

api/MWOperators

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

docs/programmers_manual/MWTree.rst

@@ -16,4 +16,3 @@ algorithm/equation/structure reasoning for having this class or where it fits wi
    :members:
    :protected-members:
    :private-members:
-

docs/programmers_manual/MultiResolutionAnalysis.rst

@@ -0,0 +1,16 @@
+------------------
+MultiResolutionAnalysis
+------------------
+
+The MultiResolutionAnalysis (MRA) class contains the methods to project objects onto the spatial grid.
+That is, to combine different functions and operators in mathematical operations, they need to be compatible; 
+they must be defined on the same computational domain and constructed using the same polynomial basis 
+(order and type). This information constitutes an MRA, which needs to be defined and passed as argument to 
+all function and operator constructors, and only functions and operators with compatible MRAs can be 
+combined in subsequent calculations.
+
+.. doxygenclass:: mrcpp::MultiResolutionAnalysis
+   :members:
+   :protected-members:
+   :private-members:
+

docs/programmers_manual/overview.rst

@@ -13,19 +13,20 @@ TODO: maybe add some low level theory/figures/algorithms before showing classes,
    OperatorNode
    OperatorTree
    BoundingBox
+   MultiResolutionAnalysis
    CrossCorrelationCache
    LegendreBasis
    InterpolatingBasis
 
    ConvolutionCalculator
 
    CornerOperatorTree
+
    TimeEvolutionOperator
    JpowerIntegrals
    special_functions
    SchrodingerEvolution_CrossCorrelation
    TimeEvolution_CrossCorrelationCalculator
    complex_apply
-
    HeatOperator
    HeatKernel

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/operators/IdentityConvolution.cpp

@@ -29,7 +29,7 @@
 
 namespace mrcpp {
 
-/** @brief Constructor  of the IdentityConvolution object
+/** @brief Constructor of the IdentityConvolution object
  *  @returns New IdentityConvolution object
  *  @param[in] mra: Which MRA the operator is defined
  *  @param[in] prec: Build precision, closeness to delta function

src/treebuilders/DerivativeCalculator.cpp

@@ -86,6 +86,34 @@ template <int D> void DerivativeCalculator<D>::printTimers() const {
     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 @@ template <int D> void DerivativeCalculator<D>::calcNode(MWNode<D> &gNode) {
     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 @@ MWNodeVector<D> DerivativeCalculator<D>::makeOperBand(const MWNode<D> &gNode, st
     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 @@ template <int D> void DerivativeCalculator<D>::applyOperator(OperatorState<D> &o
     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 @@ template <int D> double FunctionNode<D>::integrateInterpolating() const {
     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 @@ template <int D> double FunctionNode<D>::integrateInterpolating() const {
     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 @@ template <int D> void FunctionNode<D>::dealloc() {
     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 @@ template <int D>
 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 @@ template <int D> double FunctionTree<D>::integrate() const {
     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;