Tweaks to integration routine, simplier now

examples/KernelV0.cpp

@@ -33,8 +33,8 @@ int main() {
         earth->SetMagneticFieldIncDecMag( 67, 0, 52750, NANOTESLA );
 
     // Transmitter loops
-    auto Tx1 = CircularLoop(31, 15, 100, 100);
-    auto Tx2 = CircularLoop(31, 15, 100, 120);
+    auto Tx1 = CircularLoop(21, 15, 100, 100);
+    auto Tx2 = CircularLoop(21, 15, 100, 120);
     //auto Tx1 = CircularLoop(60, 15, 0, 0); // was 60
 
     auto Kern = KernelV0::NewSP();
@@ -43,13 +43,14 @@ int main() {
         Kern->SetLayeredEarthEM( earth );
         // std::cout << *Kern << std::endl;
 
-        Kern->SetIntegrationSize( (Vector3r() << 200,200,2).finished() );
-        Kern->SetIntegrationOrigin( (Vector3r() << 0,0,15).finished() );
+        Kern->SetIntegrationSize( (Vector3r() << 200,200,200).finished() );
+        Kern->SetIntegrationOrigin( (Vector3r() << 0,0,0).finished() );
         Kern->SetTolerance( 1e-13 );
+        //Kern->SetTolerance( .55 ) ; // 1%
 
         Kern->SetPulseDuration(0.020);
-        Kern->SetPulseCurrent( VectorXr::LinSpaced( 20, .01, 200 )  ); // nbins, low, high
-        Kern->SetDepthLayerInterfaces( VectorXr::LinSpaced( 20, .5, 50 ) );
+        Kern->SetPulseCurrent( VectorXr::LinSpaced( 1, 10, 200 )  ); // nbins, low, high
+        Kern->SetDepthLayerInterfaces( VectorXr::LinSpaced( 2, 5, 5.5 ) );
 
     // We could, I suppose, take the earth model in here? For non-linear that
     // may be more natural to work with?

include/KernelV0.h

@@ -265,8 +265,8 @@ namespace Lemma {
         // ====================  DATA MEMBERS  =========================
 
         int                                       nleaves;
-        int                                       minLevel=3;
-        int                                       maxLevel=10;
+        int                                       minLevel=4;
+        int                                       maxLevel=8;
 
         Real                                      VOLSUM;
         Real                                      tol=1e-11;

src/KernelV0.cpp

@@ -138,9 +138,9 @@ namespace Lemma {
         }
 
         std::cout << "Calculating K0 kernel\n";
-        MatrixXcr Kern = MatrixXcr::Zero( Interfaces.size() - 1, PulseI.size() );
-        for (int ilay=0; ilay< Interfaces.size()-1; ++ilay) {
-            for (int iq=0; iq< PulseI.size()-1; ++iq) {
+        MatrixXcr Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size() );
+        for (int ilay=0; ilay<Interfaces.size()-1; ++ilay) {
+            for (int iq=0; iq< PulseI.size(); ++iq) {
                 std::cout << "Layer " << ilay << " q " << iq << std::endl;
                 Size(2) = Interfaces(ilay+1) - Interfaces(ilay);
                 Origin(2) = Interfaces(ilay);
@@ -162,7 +162,7 @@ namespace Lemma {
     //--------------------------------------------------------------------------------------
     Complex KernelV0::IntegrateOnOctreeGrid( const int& ilay, const int& iq, bool vtkOutput) {
 
-        Vector3r cpos = (Size-Origin).array() / 2.;
+        Vector3r cpos = Origin + Size/2.;
 
         SUM = 0;
         VOLSUM = 0;
@@ -200,12 +200,15 @@ namespace Lemma {
                 kerr->SetNumberOfComponents(1);
                 kerr->SetName("nleaf");
 
+            //Real LeafVol(0);
             for (auto leaf : LeafDict) {
                 kr->InsertTuple1( leaf.first, std::real(leaf.second) );
                 ki->InsertTuple1( leaf.first, std::imag(leaf.second) );
                 km->InsertTuple1( leaf.first, std::abs(leaf.second) );
                 kid->InsertTuple1( leaf.first, leaf.first );
+                //LeafVol += std::real(leaf.second);
             }
+            //std::cout << "\n\nLeafVol=" << LeafVol << std::endl;
 
             for (auto leaf : LeafDictIdx) {
                 kerr->InsertTuple1( leaf.first, leaf.second );
@@ -427,9 +430,10 @@ namespace Lemma {
         // Next level step, interested in one level below
         // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
         Vector3r step  = size.array() / (Real)(1 << (level+1) );
-        Vector3r step2 = size.array() / (Real)(1 << (level+2) );
+        Real vol = (step(0)*step(1)*step(2));         // volume of each child
 
-        Real vol = (step2(0)*step2(1)*step2(2));     // volume of each child
+        Vector3r step2 = size.array() / (Real)(1 << (level+2) );
+        Real vol2 = (step2(0)*step2(1)*step2(2));     // volume of each child
 
         Vector3r pos =  cpos - step/2.;
         Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
@@ -479,31 +483,44 @@ namespace Lemma {
 
         Complex ksum = kvals.sum();     // Kernel sum
         // Evaluate whether or not furthur splitting is needed
-        Real err = std::abs(ksum - parentVal);
-        if ( std::abs(ksum - parentVal) > tol || level < minLevel && level < maxLevel ) {
+        if ( std::abs(ksum-parentVal) > tol ) { // || level < minLevel && level < maxLevel ) {
+        //if ( std::abs(ksum.real()-parentVal.real())>tol || std::abs(ksum.imag()-parentVal.imag()) >tol ) {
+        //if ( std::abs(ksum) > tol && level < maxLevel ) {
             oct->SubdivideLeaf(curse);
             for (int ichild=0; ichild<8; ++ichild) {
                 curse->ToChild(ichild);
                 Vector3r cp = pos; // Eigen complains about combining these
                 cp += posadd.row(ichild);
-                /* Test for position via alternative means
+                /* Test for position via alternative means */
+                /*
                 Real p[3];
                 GetPosition(curse, p);
-                std::cout << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
-                          << "\t" << cp[2] << "\t" << p[2] << "\t" <<  vol<< std::endl;
+                if ( (Vector3r(p) - cp).norm() > 1e-8 ) {
+                    std::cout << "ERROR @ nleaves" << nleaves << "\n" << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
+                              << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
+                    throw std::runtime_error("doom");
+                }
                 */
+                /* End of position test */
                 bool isleaf = EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
+                /*
                 if (isleaf) {  // Include result in final integral
-                    LeafDict[curse->GetLeafId()] = kvals(ichild);       // VTK
-                    LeafDictIdx[curse->GetLeafId()] = nleaves;          // VTK
+                    LeafDict[curse->GetLeafId()] = ksum/(8.*vol);     //kvals(ichild) / vol;     // VTK
+                    LeafDictIdx[curse->GetLeafId()] = nleaves;        // VTK
                     SUM += ksum;
                     VOLSUM += 8*vol;
                     nleaves += 1;
                 }
+                */
                 curse->ToParent();
             }
             return false;  // not leaf
         }
+        LeafDict[curse->GetLeafId()] = ksum/(8.*vol);     //kvals(ichild) / vol;     // VTK
+        LeafDictIdx[curse->GetLeafId()] = nleaves;        // VTK
+        SUM += ksum;
+        VOLSUM += 8*vol;
+        nleaves += 1;
         return true;       // leaf
     }