Lemma is an Electromagnetics API
Вы не можете выбрать более 25 тем Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251
  1. /* This file is part of Lemma, a geophysical modelling and inversion API.
  2. * More information is available at http://lemmasoftware.org
  3. */
  4. /* This Source Code Form is subject to the terms of the Mozilla Public
  5. * License, v. 2.0. If a copy of the MPL was not distributed with this
  6. * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  7. */
  8. /**
  9. * @file
  10. * @date 02/01/19 22:34:11
  11. * @version $Id$
  12. * @author Trevor Irons (ti)
  13. * @email Trevor.Irons@utah.edu
  14. * @copyright Copyright (c) 2019, University of Utah
  15. * @copyright Copyright (c) 2019, Lemma Software, LLC
  16. */
  17. #include <cxxtest/TestSuite.h>
  18. #include <FDEM1D>
  19. #include <random>
  20. #include "timer.h"
  21. #define EPSILON 1e-10
  22. using namespace Lemma;
  23. class MyAlgorithmTest : public CxxTest::TestSuite {
  24. // Put in your usual, quick unit tests here
  25. };
  26. class MyAlgorithmTestPerformance : public CxxTest::TestSuite {
  27. public:
  28. std::shared_ptr< DipoleSource > dipole;
  29. std::shared_ptr< LayeredEarthEM > earth;
  30. std::shared_ptr< FieldPoints > receivers;
  31. std::shared_ptr< EMEarth1D > EmEarth;
  32. jsw_timer timer;
  33. Real Delta;
  34. void setUp() {
  35. // Put in any code needed to set up your test,
  36. // but that should NOT be counted in the execution time.
  37. dipole = DipoleSource::NewSP();
  38. //dipole->SetType(GROUNDEDELECTRICDIPOLE);
  39. //dipole->SetPolarisation(XPOLARISATION);
  40. dipole->SetNumberOfFrequencies(1);
  41. dipole->SetMoment(1);
  42. dipole->SetFrequency(0, 4400.1000);
  43. //dipole->SetPhase(0);
  44. //dipole->SetLocation( (VectorXr(3) << 49, 49, -1e-4).finished() );
  45. dipole->SetLocation( 0, 0, -1e-4 );
  46. // Define model
  47. VectorXcr sigma(8);
  48. sigma << 0., 1e-2, .1, .01, .001, .1, .05, .2;
  49. VectorXr thick(6);
  50. thick << 10, 10, 10, 10, 10;
  51. earth = LayeredEarthEM::NewSP();
  52. earth->SetNumberOfLayers(8);
  53. earth->SetLayerConductivity(sigma);
  54. earth->SetLayerThickness(thick);
  55. receivers = FieldPoints::NewSP();
  56. Vector3r loc;
  57. Real ox = 2.;
  58. Real oy = 2.;
  59. Real oz = 2;
  60. Real dx = 0.62;
  61. Real dy = 0.62;
  62. Real dz = 0.62;
  63. int nx = 20;
  64. int ny = 20;
  65. int nz = 20;
  66. Delta = nx*ny*nz*1e-10;
  67. receivers->SetNumberOfPoints(nx*ny*nz);
  68. int ir = 0;
  69. for (int ix=0; ix<nx; ++ix) {
  70. for (int iy=0; iy<ny; ++iy) {
  71. for (int iz=0; iz<nz; ++iz) {
  72. loc << ox+ix*dx, oy+iy*dy, oz+iz*dz;
  73. //std::cout << "Receiver location " << loc.transpose() << std::endl;
  74. receivers->SetLocation(ir, loc);
  75. //oz += dz;
  76. ++ ir;
  77. }
  78. }
  79. }
  80. EmEarth = EMEarth1D::NewSP();
  81. EmEarth->SetHankelTransformMethod(CHAVE);
  82. EmEarth->SetFieldsToCalculate(BOTH); // Fortran needs this
  83. EmEarth->AttachDipoleSource(dipole);
  84. EmEarth->AttachLayeredEarthEM(earth);
  85. EmEarth->AttachFieldPoints(receivers);
  86. }
  87. void tearDown() {
  88. // Clean-up code, also NOT counted in execution time.
  89. }
  90. void test_Hz() {
  91. dipole->SetType(MAGNETICDIPOLE);
  92. dipole->SetPolarisation(ZPOLARISATION);
  93. // Put in a unit test that will be slow.
  94. std::cout << "C++\n";
  95. timer.begin();
  96. EmEarth->MakeCalc3();
  97. Real lemmaTime = timer.end();
  98. auto lc = receivers->GetEfield( 0 );
  99. #ifdef KIHALEE_EM1D
  100. receivers->ClearFields();
  101. std::cout << "\nFORTRAN KiHa\n";
  102. timer.begin();
  103. EmEarth->MakeCalc();
  104. Real kihaTime = timer.end();
  105. auto fc = receivers->GetEfield( 0 ); //0,0);
  106. std::cout << "Lemma time:" << lemmaTime << "\tKiHa time:" << kihaTime << std::endl;
  107. std::cout.precision(16);
  108. std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
  109. std::cout << "KiHa norm |" << (fc).norm() << "|" << std::endl;
  110. std::cout << "Difference norm |" << (lc - fc).norm() << "|" << std::endl;
  111. TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
  112. #endif
  113. }
  114. void test_Hx() {
  115. dipole->SetType(MAGNETICDIPOLE);
  116. dipole->SetPolarisation(XPOLARISATION);
  117. // Put in a unit test that will be slow.
  118. std::cout << "C++\n";
  119. timer.begin();
  120. EmEarth->MakeCalc3();
  121. Real lemmaTime = timer.end();
  122. auto lc = receivers->GetEfield( 0 );
  123. #ifdef KIHALEE_EM1D
  124. receivers->ClearFields();
  125. std::cout << "\nFORTRAN KiHa\n";
  126. timer.begin();
  127. EmEarth->MakeCalc();
  128. Real kihaTime = timer.end();
  129. auto fc = receivers->GetEfield( 0 ); //0,0);
  130. std::cout << "Lemma time:" << lemmaTime << "\tKiHa time:" << kihaTime << std::endl;
  131. std::cout.precision(16);
  132. std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
  133. std::cout << "KiHa norm |" << (fc).norm() << "|" << std::endl;
  134. std::cout << "Difference norm |" << (lc - fc).norm() << "|" << std::endl;
  135. TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
  136. #endif
  137. }
  138. void test_Hy() {
  139. dipole->SetType(MAGNETICDIPOLE);
  140. dipole->SetPolarisation(YPOLARISATION);
  141. // Put in a unit test that will be slow.
  142. std::cout << "C++\n";
  143. timer.begin();
  144. EmEarth->MakeCalc3();
  145. Real lemmaTime = timer.end();
  146. auto lc = receivers->GetEfield( 0 );
  147. #ifdef KIHALEE_EM1D
  148. receivers->ClearFields();
  149. std::cout << "\nFORTRAN KiHa\n";
  150. timer.begin();
  151. EmEarth->MakeCalc();
  152. Real kihaTime = timer.end();
  153. auto fc = receivers->GetEfield( 0 ); //0,0);
  154. std::cout << "Lemma time:" << lemmaTime << "\tKiHa time:" << kihaTime << std::endl;
  155. std::cout.precision(16);
  156. std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
  157. std::cout << "KiHa norm |" << (fc).norm() << "|" << std::endl;
  158. std::cout << "Difference norm |" << (lc - fc).norm() << "|" << std::endl;
  159. TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
  160. #endif
  161. }
  162. void test_Ex() {
  163. dipole->SetType(GROUNDEDELECTRICDIPOLE);
  164. dipole->SetPolarisation(XPOLARISATION);
  165. // Put in a unit test that will be slow.
  166. std::cout << "C++\n";
  167. timer.begin();
  168. EmEarth->MakeCalc3();
  169. Real lemmaTime = timer.end();
  170. auto lc = receivers->GetEfield( 0 );
  171. #ifdef KIHALEE_EM1D
  172. receivers->ClearFields();
  173. std::cout << "\nFORTRAN KiHa\n";
  174. timer.begin();
  175. EmEarth->MakeCalc();
  176. Real kihaTime = timer.end();
  177. auto fc = receivers->GetEfield( 0 ); //0,0);
  178. std::cout << "Lemma time:" << lemmaTime << "\tKiHa time:" << kihaTime << std::endl;
  179. std::cout.precision(16);
  180. std::cout << "Lemma norm |" << (lc).norm() << "|" << std::endl;
  181. std::cout << "KiHa norm |" << (fc).norm() << "|" << std::endl;
  182. std::cout << "Difference norm |" << (lc - fc).norm() << "|" << std::endl;
  183. TS_ASSERT_DELTA((lc-fc).norm(), 0.0, Delta);
  184. #endif
  185. }
  186. };