Clean up of KiHa test

Modules/FDEM1D/examples/EMDipEarth1D.cpp

@@ -68,7 +68,7 @@ int main() {
 		Real oy =    20.;
 		Real depth = 18.10;
 		Real dz = 2.6;
-		int  nz = 1;
+		int  nz = 10000;
 
 		receivers->SetNumberOfPoints(nz);
 		int ir = 0;
@@ -83,6 +83,8 @@ int main() {
 
     auto EmEarth = EMEarth1D::NewSP();
         //EmEarth->SetHankelTransformMethod(DIGITALFILTERING);
+        EmEarth->SetHankelTransformMethod(CHAVE);
+
         EmEarth->SetFieldsToCalculate(BOTH); // Fortran needs this
 		EmEarth->AttachDipoleSource(dipole);
 		EmEarth->AttachLayeredEarthEM(earth);

Modules/FDEM1D/examples/inp/config.inp

@@ -1,13 +1,13 @@
 // Hankel Transform type uncomment desired
-FHTKEY51       // Key's   51 point FHT
-//FHTKEY101      // Key's  101 point FHT
-//FHTKEY201      // Key's  201 point FHT
+CHAVE       // Key's   51 point FHT
+//CHAVE      // Key's  101 point FHT
+//CHAVE      // Key's  201 point FHT
 //FHTKONG61      // Kong's  61 point FHT
 //FHTKONG121     // Kong's 121 point FHT
 //FHTKONG241	 // Kong's 241 point FHT
-//ANDERSON801    // Anderson's 801 Point FHT, with lagged convolution
+//CHAVE    // Anderson's 801 Point FHT, with lagged convolution
 //CHAVE          // Chave's Gaussian Quadrature
-//QWEKEY         // Key's Gaussian Qwuadrature with extraploation
+//CHAVE         // Key's Gaussian Qwuadrature with extraploation
 //IRONS
 .1               // minimum dipole ratio
 1e-5             // minumum dipole moment

Modules/FDEM1D/examples/plottimings.py

@@ -56,7 +56,7 @@ for compiler in BENCH:
 
 plt.gca().set_yscale('log')    
 plt.gca().set_title( str(compiler) + " " + str(BENCH[compiler]["version"]) ) 
-plt.suptitle( cpuinfo.get_cpu_info()['brand'] ) 
+#plt.suptitle( cpuinfo.get_cpu_info()['brand'] ) 
 
 
 

Modules/FDEM1D/src/EMEarth1D.cpp

@@ -680,7 +680,7 @@ namespace Lemma {
         Real rho = (Receivers->GetLocation(irec).head<2>() - tDipole->GetLocation().head<2>()).norm();
         //Real rho = ( ((Receivers->GetLocation(irec) - tDipole->GetLocation()).head(2)).eval() ).norm();
 
-        tDipole->SetKernels(ifreq, FieldsToCalculate, Receivers, irec, Earth);
+        tDipole->SetKernels( ifreq, FieldsToCalculate, Receivers, irec, Earth );
         Hankel->ComputeRelated( rho, tDipole->GetKernelManager() );
         tDipole->UpdateFields( ifreq,  Hankel, wavef );
     }
@@ -755,6 +755,7 @@ namespace Lemma {
         #pragma omp parallel
         #endif
         { // OpenMP Parallel Block
+
             #ifdef LEMMAUSEOMP
             int tid = omp_get_thread_num();
             int nthreads = omp_get_num_threads();
@@ -763,6 +764,7 @@ namespace Lemma {
             int nthreads=1;
             #endif
             auto tDipole = Dipole->Clone();
+
             std::shared_ptr<HankelTransform> Hankel;
             switch (HankelType) {
                 case ANDERSON801:
@@ -787,6 +789,7 @@ namespace Lemma {
                     std::cerr << "Hankel transform cannot be created\n";
                     exit(EXIT_FAILURE);
             }
+
             if ( tDipole->GetNumberOfFrequencies() < Receivers->GetNumberOfPoints() ) {
                 for (int ifreq=0; ifreq<tDipole->GetNumberOfFrequencies(); ++ifreq) {
                     // Propogation constant in free space being input to Hankel

Modules/FDEM1D/src/FHTAnderson801.cpp

@@ -1042,7 +1042,7 @@ namespace Lemma {
 		//Zans.resize(this->NumConv, (int)(this->kernelVec.size()));
 		//Zans.setZero();
         //Zwork = Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(801, (int)(this->kernelVec.size()));
-        Zans  = Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(this->NumConv, (int)(this->kernelVec.size()));
+        Zans = Eigen::Matrix<Complex, Eigen::Dynamic, Eigen::Dynamic>::Zero(this->NumConv, (int)(this->kernelVec.size()));
 
  		// 1010 Loop
  		for (int ilag=0; ilag < this->NumConv; ++ilag) {

Modules/FDEM1D/src/KernelEM1DSpec.cpp

@@ -1682,7 +1682,6 @@ namespace Lemma {
     }
 
     // TODO in PotentialBelowSourceLayer:
-
     template<>
     Complex KernelEM1DSpec<TM, 0, INAIR, INGROUND>::PotentialBelowSourceLayer(const Real &ra) {
         Complex dd =  ((Real)(1.)+ReflCalc->rtd(1)*ReflCalc->cf(1));
@@ -1815,7 +1814,7 @@ namespace Lemma {
             }
             p += (ReflCalc->u(n)-ut) * ReflCalc->LayerDepth(n-1);
         }
-        Complex con = SR_SN(0, 0) * std::exp(ReflCalc->uk*ReflCalc->tx_z - ReflCalc->um*ReflCalc->rx_z+ p);
+        Complex con = SR_SN(0, 0) * std::exp(ReflCalc->uk*ReflCalc->tx_z - ReflCalc->um*ReflCalc->rx_z + p);
         if (ReflCalc->layr < ReflCalc->Earth->GetNumberOfLayers()-1) {
             con += SR_SN(0, 2) * ReflCalc->rtd(ReflCalc->layr) * std::exp(ReflCalc->uk*ReflCalc->tx_z-
                         ReflCalc->um*((Real)(2.)*ReflCalc->LayerDepth(ReflCalc->layr)-ReflCalc->rx_z)+p);

Modules/FDEM1D/testing/BenchKiHa.h

@@ -80,9 +80,9 @@ public:
 		Real dy = 20;
 		Real dz = 20;
 
-		int  nx = 5;  //13
-		int  ny = 5;  //13
-		int  nz = 5;  //13
+		int  nx = 11;  //13
+		int  ny = 11;  //13
+		int  nz = 11;  //13
         Delta = nx*ny*nz*1e-9;
 
 		receivers->SetNumberOfPoints(nx*ny*nz);
@@ -114,40 +114,41 @@ public:
 
    void test_Hz() {
 
+        std::cout.precision(4);
+
         dipole->SetType(MAGNETICDIPOLE);
 		dipole->SetPolarisation(ZPOLARISATION);
 
         // Put in a unit test that will be slow.
         std::cout << "MAGNETICDIPOLE Z polarisation" << std::endl;
         std::cout << "=====================================\n";
-	    std::cout << "Lemma/C++: ";
+	    std::cout << std::setw(18) << "Lemma/C++: ";
         std::cout.flush();
 
   	    timer.begin();
 	    EmEarth->MakeCalc3();
 	    Real lemmaTime = timer.end();
-        std::cout << lemmaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << lemmaTime << std::setw(6) << " [s]" << std::endl;
 
         auto lc = receivers->GetEfield( 0 );
 
         #ifdef KIHALEE_EM1D
 	    receivers->ClearFields();
-        std::cout << "FORTRAN KiHa: ";
+        std::cout << std::setw(18) << "KiHa/Fortran: ";
         std::cout.flush();
 
   	    timer.begin();
  	    EmEarth->MakeCalc();
 	    Real kihaTime = timer.end();
-
-        std::cout << kihaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << kihaTime << std::setw(6) << " [s]" << std::endl;
 
         auto fc = receivers->GetEfield( 0 ); //0,0);
 
-        std::cout.precision(16);
-        std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
-        std::cout << "KiHa norm  |" << (fc).norm() << "|" << std::endl;
-        std::cout << "Difference norm |" << (lc - fc).norm() << "|\n";
-        std::cout << "Speedup:" << kihaTime/lemmaTime << "\n" << std::endl;
+        //std::cout.precision(16);
+        std::cout << std::setw(18) << "Lemma norm: "      << std::setw(14) << (lc).norm() << std::endl;
+        std::cout << std::setw(18) << "KiHa norm: "       << std::setw(14) <<  (fc).norm() << std::endl;
+        std::cout << std::setw(18) << "Difference norm: " << std::setw(14) << (lc - fc).norm() << "\n";
+        std::cout << std::setw(18) << "Speedup: "         << std::setw(14) << kihaTime/lemmaTime << "\n" << std::endl;
 
         TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
         #endif
@@ -155,37 +156,39 @@ public:
 
    void test_Hx() {
 
+        std::cout.precision(4);
+
         dipole->SetType(MAGNETICDIPOLE);
 		dipole->SetPolarisation(XPOLARISATION);
 
         // Put in a unit test that will be slow.
         std::cout << "MAGNETICDIPOLE X polarisation" << std::endl;
         std::cout << "=====================================\n";
-	    std::cout << "Lemma/C++: ";
+	    std::cout << std::setw(18) << "Lemma/C++: ";
         std::cout.flush();
 
   	    timer.begin();
 	    EmEarth->MakeCalc3();
 	    Real lemmaTime = timer.end();
-        std::cout << lemmaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << lemmaTime << std::setw(6) << " [s]" << std::endl;
 
         auto lc = receivers->GetEfield( 0 );
 
         #ifdef KIHALEE_EM1D
 	    receivers->ClearFields();
-        std::cout << "FORTRAN KiHa: ";
+        std::cout << std::setw(18) << "KiHa/Fortran: ";
   	    timer.begin();
  	    EmEarth->MakeCalc();
 	    Real kihaTime = timer.end();
-        std::cout << kihaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << kihaTime << std::setw(6) << " [s]" << std::endl;
 
         auto fc = receivers->GetEfield( 0 ); //0,0);
 
-        std::cout.precision(16);
-        std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
-        std::cout << "KiHa norm  |" << (fc).norm() << "|" << std::endl;
-        std::cout << "Difference norm |" << (lc - fc).norm() << "|\n";
-        std::cout << "Speedup:" << kihaTime/lemmaTime << "\n" << std::endl;
+        //std::cout.precision(16);
+        std::cout << std::setw(18) << "Lemma norm: "      << std::setw(14) << (lc).norm() << std::endl;
+        std::cout << std::setw(18) << "KiHa norm: "       << std::setw(14) <<  (fc).norm() << std::endl;
+        std::cout << std::setw(18) << "Difference norm: " << std::setw(14) << (lc - fc).norm() << "\n";
+        std::cout << std::setw(18) << "Speedup: "         << std::setw(14) << kihaTime/lemmaTime << "\n" << std::endl;
 
         TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
         #endif
@@ -193,38 +196,40 @@ public:
 
    void test_Hy() {
 
+        std::cout.precision(4);
+
         dipole->SetType(MAGNETICDIPOLE);
 		dipole->SetPolarisation(YPOLARISATION);
 
         // Put in a unit test that will be slow.
         std::cout << "MAGNETICDIPOLE Y polarisation" << std::endl;
         std::cout << "=====================================\n";
-	    std::cout << "Lemma/C++: ";
+	    std::cout << std::setw(18) << "Lemma/C++: ";
         std::cout.flush();
 
   	    timer.begin();
 	    EmEarth->MakeCalc3();
 	    Real lemmaTime = timer.end();
-        std::cout << lemmaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << lemmaTime << std::setw(6) << " [s]" << std::endl;
 
         auto lc = receivers->GetEfield( 0 );
 
         #ifdef KIHALEE_EM1D
 	    receivers->ClearFields();
-        std::cout << "FORTRAN KiHa: ";
+        std::cout << std::setw(18) << "KiHa/Fortran: ";
         std::cout.flush();
   	    timer.begin();
  	    EmEarth->MakeCalc();
 	    Real kihaTime = timer.end();
-        std::cout << kihaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << kihaTime << std::setw(6) << " [s]" << std::endl;
 
         auto fc = receivers->GetEfield( 0 ); //0,0);
 
-        std::cout.precision(16);
-        std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
-        std::cout << "KiHa norm  |" << (fc).norm() << "|" << std::endl;
-        std::cout << "Difference norm |" << (lc - fc).norm() << "|\n";
-        std::cout << "Speedup:" << kihaTime/lemmaTime << "\n" << std::endl;
+        //std::cout.precision(16);
+        std::cout << std::setw(18) << "Lemma norm: "      << std::setw(14) << (lc).norm() << std::endl;
+        std::cout << std::setw(18) << "KiHa norm: "       << std::setw(14) <<  (fc).norm() << std::endl;
+        std::cout << std::setw(18) << "Difference norm: " << std::setw(14) << (lc - fc).norm() << "\n";
+        std::cout << std::setw(18) << "Speedup: "         << std::setw(14) << kihaTime/lemmaTime << "\n" << std::endl;
 
         TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
         #endif
@@ -232,38 +237,40 @@ public:
 
    void test_Ex() {
 
+        std::cout.precision(4);
+
         dipole->SetType(GROUNDEDELECTRICDIPOLE);
 		dipole->SetPolarisation(XPOLARISATION);
 
         // Put in a unit test that will be slow.
         std::cout << "GROUNDEDELECTRICDIPOLE X polarisation" << std::endl;
         std::cout << "=====================================\n";
-	    std::cout << "Lemma/C++: ";
+	    std::cout << std::setw(18) << "Lemma/C++: ";
         std::cout.flush();
 
   	    timer.begin();
 	    EmEarth->MakeCalc3();
 	    Real lemmaTime = timer.end();
-        std::cout << lemmaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << lemmaTime << std::setw(6) << " [s]" << std::endl;
 
         auto lc = receivers->GetEfield( 0 );
 
         #ifdef KIHALEE_EM1D
 	    receivers->ClearFields();
-        std::cout << "FORTRAN KiHa: ";
+        std::cout << std::setw(18) << "KiHa/Fortran: ";
         std::cout.flush();
   	    timer.begin();
  	    EmEarth->MakeCalc();
 	    Real kihaTime = timer.end();
-        std::cout << kihaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << kihaTime << std::setw(6) << " [s]" << std::endl;
 
         auto fc = receivers->GetEfield( 0 ); //0,0);
 
-        std::cout.precision(16);
-        std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
-        std::cout << "KiHa norm  |" << (fc).norm() << "|" << std::endl;
-        std::cout << "Difference norm |" << (lc - fc).norm() << "|\n";
-        std::cout << "Speedup:" << kihaTime/lemmaTime << "\n" << std::endl;
+        //std::cout.precision(16);
+        std::cout << std::setw(18) << "Lemma norm: "      << std::setw(14) << (lc).norm() << std::endl;
+        std::cout << std::setw(18) << "KiHa norm: "       << std::setw(14) <<  (fc).norm() << std::endl;
+        std::cout << std::setw(18) << "Difference norm: " << std::setw(14) << (lc - fc).norm() << "\n";
+        std::cout << std::setw(18) << "Speedup: "         << std::setw(14) << kihaTime/lemmaTime << "\n" << std::endl;
 
         TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
         #endif
@@ -271,41 +278,41 @@ public:
 
    void test_Ey() {
 
+        std::cout.precision(4);
+
         dipole->SetType(GROUNDEDELECTRICDIPOLE);
 		dipole->SetPolarisation(YPOLARISATION);
 
         // Put in a unit test that will be slow.
         std::cout << "GROUNDEDELECTRICDIPOLE Y polarisation" << std::endl;
         std::cout << "=====================================\n";
-	    std::cout << "Lemma/C++: ";
+	    std::cout << std::setw(18) << "Lemma/C++: ";
         std::cout.flush();
 
   	    timer.begin();
 	    EmEarth->MakeCalc3();
 	    Real lemmaTime = timer.end();
-        std::cout << lemmaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << lemmaTime << std::setw(6) << " [s]" << std::endl;
 
         auto lc = receivers->GetEfield( 0 );
 
         #ifdef KIHALEE_EM1D
 	    receivers->ClearFields();
-        std::cout << "KiHa/FORTRAN: ";
+        std::cout << std::setw(18) << "KiHa/Fortran: ";
         std::cout.flush();
 
         timer.begin();
  	    EmEarth->MakeCalc();
 	    Real kihaTime = timer.end();
-        std::cout << kihaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << kihaTime << std::setw(6) << " [s]" << std::endl;
 
         auto fc = receivers->GetEfield( 0 ); //0,0);
 
-        std::cout << "Lemma time:" << lemmaTime << "\tKiHa time:" << kihaTime << std::endl;
-
-        std::cout.precision(16);
-        std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
-        std::cout << "KiHa norm  |" << (fc).norm() << "|" << std::endl;
-        std::cout << "Difference norm |" << (lc - fc).norm() << "|\n";
-        std::cout << "Speedup:" << kihaTime/lemmaTime << "\n" << std::endl;
+        //std::cout.precision(16);
+        std::cout << std::setw(18) << "Lemma norm: "      << std::setw(14) << (lc).norm() << std::endl;
+        std::cout << std::setw(18) << "KiHa norm: "       << std::setw(14) <<  (fc).norm() << std::endl;
+        std::cout << std::setw(18) << "Difference norm: " << std::setw(14) << (lc - fc).norm() << "\n";
+        std::cout << std::setw(18) << "Speedup: "         << std::setw(14) << kihaTime/lemmaTime << "\n" << std::endl;
 
         TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
         #endif
@@ -313,41 +320,41 @@ public:
 
    void test_Ez() {
 
+        std::cout.precision(4);
+
         dipole->SetType(GROUNDEDELECTRICDIPOLE);
 		dipole->SetPolarisation(ZPOLARISATION);
 
         // Put in a unit test that will be slow.
         std::cout << "GROUNDEDELECTRICDIPOLE Z polarisation" << std::endl;
         std::cout << "=====================================\n";
-	    std::cout << "Lemma/C++: ";
+	    std::cout << std::setw(18) << "Lemma/C++: ";
         std::cout.flush();
 
   	    timer.begin();
 	    EmEarth->MakeCalc3();
 	    Real lemmaTime = timer.end();
-        std::cout << lemmaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << lemmaTime << std::setw(6) << " [s]" << std::endl;
 
         auto lc = receivers->GetEfield( 0 );
 
         #ifdef KIHALEE_EM1D
 	    receivers->ClearFields();
-        std::cout << "KiHa/FORTRAN: ";
+        std::cout << std::setw(18) << "KiHa/Fortran: ";
         std::cout.flush();
 
         timer.begin();
  	    EmEarth->MakeCalc();
 	    Real kihaTime = timer.end();
-        std::cout << kihaTime << " [s]" << std::endl;
+        std::cout << std::setw(14) << kihaTime << std::setw(6) << " [s]" << std::endl;
 
         auto fc = receivers->GetEfield( 0 ); //0,0);
 
-        std::cout << "Lemma time:" << lemmaTime << "\tKiHa time:" << kihaTime << std::endl;
-
-        std::cout.precision(16);
-        std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
-        std::cout << "KiHa norm  |" << (fc).norm() << "|" << std::endl;
-        std::cout << "Difference norm |" << (lc - fc).norm() << "|\n";
-        std::cout << "Speedup:" << kihaTime/lemmaTime << "\n" << std::endl;
+        //std::cout.precision(16);
+        std::cout << std::setw(18) << "Lemma norm: "      << std::setw(14) << (lc).norm() << std::endl;
+        std::cout << std::setw(18) << "KiHa norm: "       << std::setw(14) <<  (fc).norm() << std::endl;
+        std::cout << std::setw(18) << "Difference norm: " << std::setw(14) << (lc - fc).norm() << "\n";
+        std::cout << std::setw(18) << "Speedup: "         << std::setw(14) << kihaTime/lemmaTime << "\n" << std::endl;
 
         TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
         #endif