KiHa unit test expanded

Modules/FDEM1D/CMakeLists.txt

@@ -15,7 +15,13 @@ target_link_libraries(fdem1d "lemmacore")
 
 if ( KIHA_EM1D )
 	target_link_libraries(fdem1d "kihaem1d")
-	target_link_libraries(fdem1d "gfortran")
+	if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+  		# using GCC
+		target_link_libraries(fdem1d "gfortran")
+	elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+  		# using Intel C++
+		target_link_libraries(fdem1d "ifcore" "imf" "svml" "intlc")
+	endif()
 endif()
 
 # Testing

Modules/FDEM1D/examples/Bench.sh

@@ -0,0 +1,50 @@
+#!/bin/bash
+
+# this is a simple benchmarking script for Hantenna 
+
+for i in {1..12}
+do
+   echo "Welcome $i times"
+   export OMP_NUM_THREADS=$i
+   ./Hantenna inp/trans.inp inp/cond.inp inp/points.inp inp/config.inp
+done
+
+sed -i -- 's/FHTKEY51/FHTKEY101/g' inp/config.inp 
+for i in {1..12}
+do
+   echo "Welcome $i times"
+   export OMP_NUM_THREADS=$i
+   ./Hantenna inp/trans.inp inp/cond.inp inp/points.inp inp/config.inp
+done
+
+sed -i -- 's/FHTKEY101/FHTKEY201/g' inp/config.inp 
+for i in {1..12}
+do
+   echo "Welcome $i times"
+   export OMP_NUM_THREADS=$i
+   ./Hantenna inp/trans.inp inp/cond.inp inp/points.inp inp/config.inp
+done
+
+sed -i -- 's/FHTKEY201/ANDERSON801/g' inp/config.inp 
+for i in {1..12}
+do
+   echo "Welcome $i times"
+   export OMP_NUM_THREADS=$i
+   ./Hantenna inp/trans.inp inp/cond.inp inp/points.inp inp/config.inp
+done
+
+sed -i -- 's/ANDERSON801/QWEKEY/g' inp/config.inp 
+for i in {1..12}
+do
+   echo "Welcome $i times"
+   export OMP_NUM_THREADS=$i
+   ./Hantenna inp/trans.inp inp/cond.inp inp/points.inp inp/config.inp
+done
+
+sed -i -- 's/QWEKEY/CHAVE/g' inp/config.inp 
+for i in {1..12}
+do
+   echo "Welcome $i times"
+   export OMP_NUM_THREADS=$i
+   ./Hantenna inp/trans.inp inp/cond.inp inp/points.inp inp/config.inp
+done

Modules/FDEM1D/examples/CMakeLists.txt

@@ -41,3 +41,4 @@ INSTALL_TARGETS( "/share/FDEM1D/"
 )
 
 install (DIRECTORY inp  DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FDEM1D/" )
+install (FILES  Bench.sh plottimings.py  DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FDEM1D/" )

Modules/FDEM1D/examples/Hantenna.cpp

@@ -213,6 +213,14 @@ const char *buildString = __DATE__ ", " __TIME__;
     }
     hreal.close();
     // Clean up
+
+    // report timings
+    #ifdef LEMMAUSEOMP
+    std::ofstream outfile;
+    outfile.open("timings.csv", std::ios_base::app);
+    outfile << compiler << "," << ver << "," <<  buildString << "," << config[0] << "," <<  omp_get_max_threads() << "," << paTime/60. << "\n";
+    #endif
+
 }
 
 std::vector<Real>  readinpfile(const std::string& fname) {

Modules/FDEM1D/examples/inp/config.inp

@@ -1,9 +1,9 @@
 // Hankel Transform type uncomment desired
-//FHTKEY51       // Key's   51 point FHT
+FHTKEY51       // Key's   51 point FHT
 //FHTKEY101      // Key's  101 point FHT
 //FHTKEY201      // Key's  201 point FHT
 //FHTKONG61      // Kong's  61 point FHT
-FHTKONG121     // Kong's 121 point FHT
+//FHTKONG121     // Kong's 121 point FHT
 //FHTKONG241	 // Kong's 241 point FHT
 //ANDERSON801    // Anderson's 801 Point FHT, with lagged convolution
 //CHAVE          // Chave's Gaussian Quadrature

Modules/FDEM1D/examples/plottimings.py

@@ -0,0 +1,69 @@
+import matplotlib.pyplot as plt 
+import numpy as np 
+from collections import defaultdict
+import cpuinfo # pip install py-cpuinfo
+
+light_grey = np.array([float(225)/float(255)]*3)
+
+def fixLeg(legend):
+    rect = legend.get_frame()
+    #rect.set_color('None')
+    rect.set_facecolor(light_grey)
+    rect.set_linewidth(0.0)
+    rect.set_alpha(1.0)
+
+def deSpine(ax1):
+    spines_to_remove = ['top', 'right']
+    for spine in spines_to_remove:
+        ax1.spines[spine].set_visible(False)
+    #ax1.xaxis.set_ticks_position('none')
+    #ax1.yaxis.set_ticks_position('none')
+    ax1.get_xaxis().tick_bottom()
+    ax1.get_yaxis().tick_left()
+
+timings = np.recfromtxt(r"timings.csv", delimiter=r",", encoding=None)
+BENCH = {} # defaultdict(str)
+
+for i,line in enumerate(timings): 
+
+    # check if compiler has been added already 
+    if line[0].strip() in BENCH:
+        pass
+    else:
+        BENCH[line[0].strip()] = {}
+        print(line[1])
+        BENCH[line[0].strip()]["version"] = line[1].strip()
+
+    # check to see if Hankel has been added already 
+    if line[4].strip() in BENCH[line[0].strip()]:
+        pass
+    else: 
+        BENCH[line[0].strip()][line[4].strip()] = {}
+        BENCH[line[0].strip()][line[4].strip()]["NP"] = []
+        BENCH[line[0].strip()][line[4].strip()]["Time"] = []
+        
+    BENCH[line[0].strip()][line[4].strip()]["NP"].append( line[5] )
+    BENCH[line[0].strip()][line[4].strip()]["Time"].append( line[6] )
+
+#print (BENCH)
+fig = plt.figure()
+
+for compiler in BENCH:    
+    print (compiler)
+    for hankel in BENCH[compiler]:   
+        if hankel != "version": 
+            plt.plot(BENCH[compiler][hankel]["NP"], BENCH[compiler][hankel]["Time"], ".-", label=str(hankel))
+
+plt.gca().set_yscale('log')    
+plt.gca().set_title( str(compiler) + " " + str(BENCH[compiler]["version"]) ) 
+plt.suptitle( cpuinfo.get_cpu_info()['brand'] ) 
+
+
+
+leg = plt.legend()
+
+deSpine(plt.gca())
+fixLeg(leg)
+
+plt.savefig("timings.png")
+plt.show()

Modules/FDEM1D/testing/BenchKiHa.h

@@ -0,0 +1,251 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+/**
+ * @file
+ * @date      02/01/19 22:34:11
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@utah.edu
+ * @copyright Copyright (c) 2019, University of Utah
+ * @copyright Copyright (c) 2019, Lemma Software, LLC
+ */
+
+#include <cxxtest/TestSuite.h>
+#include <FDEM1D>
+#include <random>
+#include "timer.h"
+
+#define EPSILON 1e-10
+using namespace Lemma;
+
+class MyAlgorithmTest : public CxxTest::TestSuite {
+   // Put in your usual, quick unit tests here
+};
+
+class MyAlgorithmTestPerformance : public CxxTest::TestSuite {
+
+public:
+
+    std::shared_ptr< DipoleSource > dipole;
+    std::shared_ptr< LayeredEarthEM > earth;
+    std::shared_ptr< FieldPoints > receivers;
+    std::shared_ptr< EMEarth1D > EmEarth;
+
+	jsw_timer timer;
+
+    Real Delta;
+
+    void setUp() {
+        // Put in any code needed to set up your test,
+        // but that should NOT be counted in the execution time.
+	    dipole = DipoleSource::NewSP();
+
+        //dipole->SetType(GROUNDEDELECTRICDIPOLE);
+		//dipole->SetPolarisation(XPOLARISATION);
+
+
+		dipole->SetNumberOfFrequencies(1);
+		dipole->SetMoment(1);
+		dipole->SetFrequency(0, 4400.1000);
+		//dipole->SetPhase(0);
+		//dipole->SetLocation( (VectorXr(3) << 49, 49, -1e-4).finished() );
+		dipole->SetLocation( 0, 0, -1e-4  );
+
+	    // Define model
+	    VectorXcr sigma(2);
+		sigma << 0., 1e-3;//, .1;//, .01, .001;
+	    VectorXr  thick(1);
+		thick << 10;//, 10, 10;
+
+	    earth = LayeredEarthEM::NewSP();
+        earth->SetNumberOfLayers(2);
+		earth->SetLayerConductivity(sigma);
+		//earth->SetLayerThickness(thick);
+
+	    receivers = FieldPoints::NewSP();
+		Vector3r loc;
+
+        Real ox = 2.;
+		Real oy = 2.;
+		Real oz = 2;
+
+		Real dx = 0.62;
+		Real dy = 0.62;
+		Real dz = 0.62;
+
+		int  nx = 20;
+		int  ny = 20;
+		int  nz = 20;
+        Delta = nx*ny*nz*1e-10;
+
+		receivers->SetNumberOfPoints(nx*ny*nz);
+		int ir = 0;
+		for (int ix=0; ix<nx; ++ix) {
+		for (int iy=0; iy<ny; ++iy) {
+		for (int iz=0; iz<nz; ++iz) {
+			loc << ox+ix*dx, oy+iy*dy, oz+iz*dz;
+            //std::cout << "Receiver location " << loc.transpose() << std::endl;
+			receivers->SetLocation(ir, loc);
+			//oz += dz;
+			++ ir;
+		}
+        }
+        }
+
+        EmEarth = EMEarth1D::NewSP();
+        EmEarth->SetHankelTransformMethod(CHAVE);
+        EmEarth->SetFieldsToCalculate(BOTH); // Fortran needs this
+		EmEarth->AttachDipoleSource(dipole);
+		EmEarth->AttachLayeredEarthEM(earth);
+		EmEarth->AttachFieldPoints(receivers);
+
+   }
+
+   void tearDown() {
+      // Clean-up code, also NOT counted in execution time.
+   }
+
+   void test_Hz() {
+
+        dipole->SetType(MAGNETICDIPOLE);
+		dipole->SetPolarisation(ZPOLARISATION);
+
+        // Put in a unit test that will be slow.
+	    std::cout << "C++\n";
+
+  	    timer.begin();
+	    EmEarth->MakeCalc3();
+	    Real lemmaTime = timer.end();
+
+        auto lc = receivers->GetEfield( 0 );
+
+        #ifdef KIHALEE_EM1D
+	    receivers->ClearFields();
+        std::cout << "\nFORTRAN KiHa\n";
+  	    timer.begin();
+ 	    EmEarth->MakeCalc();
+	    Real kihaTime = timer.end();
+
+        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() << "|" << std::endl;
+
+        TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
+        #endif
+   }
+
+   void test_Hx() {
+
+        dipole->SetType(MAGNETICDIPOLE);
+		dipole->SetPolarisation(XPOLARISATION);
+
+        // Put in a unit test that will be slow.
+	    std::cout << "C++\n";
+
+  	    timer.begin();
+	    EmEarth->MakeCalc3();
+	    Real lemmaTime = timer.end();
+
+        auto lc = receivers->GetEfield( 0 );
+
+        #ifdef KIHALEE_EM1D
+	    receivers->ClearFields();
+        std::cout << "\nFORTRAN KiHa\n";
+  	    timer.begin();
+ 	    EmEarth->MakeCalc();
+	    Real kihaTime = timer.end();
+
+        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() << "|" << std::endl;
+
+        TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
+        #endif
+   }
+
+   void test_Hy() {
+
+        dipole->SetType(MAGNETICDIPOLE);
+		dipole->SetPolarisation(YPOLARISATION);
+
+        // Put in a unit test that will be slow.
+	    std::cout << "C++\n";
+
+  	    timer.begin();
+	    EmEarth->MakeCalc3();
+	    Real lemmaTime = timer.end();
+
+        auto lc = receivers->GetEfield( 0 );
+
+        #ifdef KIHALEE_EM1D
+	    receivers->ClearFields();
+        std::cout << "\nFORTRAN KiHa\n";
+  	    timer.begin();
+ 	    EmEarth->MakeCalc();
+	    Real kihaTime = timer.end();
+
+        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() << "|" << std::endl;
+
+        TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
+        #endif
+   }
+
+   void test_Ex() {
+
+        dipole->SetType(GROUNDEDELECTRICDIPOLE);
+		dipole->SetPolarisation(XPOLARISATION);
+
+        // Put in a unit test that will be slow.
+	    std::cout << "C++\n";
+
+  	    timer.begin();
+	    EmEarth->MakeCalc3();
+	    Real lemmaTime = timer.end();
+
+        auto lc = receivers->GetEfield( 0 );
+
+        #ifdef KIHALEE_EM1D
+	    receivers->ClearFields();
+        std::cout << "\nFORTRAN KiHa\n";
+  	    timer.begin();
+ 	    EmEarth->MakeCalc();
+	    Real kihaTime = timer.end();
+
+        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() << "|" << std::endl;
+
+        TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
+        #endif
+   }
+
+};

Modules/FDEM1D/testing/CMakeLists.txt

@@ -6,9 +6,13 @@ target_link_libraries(unittest_FEM1D_GetNameCheck "lemmacore" "fdem1d" "yaml-cpp
 CXXTEST_ADD_TEST(unittest_FEM1D_SerializeCheck SerializeCheck.cc ${CMAKE_CURRENT_SOURCE_DIR}/SerializeCheck.h)
 target_link_libraries(unittest_FEM1D_SerializeCheck "lemmacore" "fdem1d" "yaml-cpp")
 
+CXXTEST_ADD_TEST(benchKiHa BenchKiHa.cc ${CMAKE_CURRENT_SOURCE_DIR}/BenchKiHa.h)
+target_link_libraries(benchKiHa "lemmacore" "fdem1d" "yaml-cpp")
+
 set_target_properties(unittest_FEM1D_GetNameCheck
-	              unittest_FEM1D_SerializeCheck 
-       	PROPERTIES 
+	                  unittest_FEM1D_SerializeCheck 
+                      benchKiHa
+    PROPERTIES 
 	CXX_STANDARD 14
 	CXX_STANDARD_REQUIRED ON
-	)
+)