Surface NMR forward modelling
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

Coupling.cpp 19KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451
  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 11/11/2016 01:47:25 PM
  11. * @author Trevor Irons (ti)
  12. * @email tirons@egi.utah.edu
  13. * @copyright Copyright (c) 2016, University of Utah
  14. * @copyright Copyright (c) 2016, Lemma Software, LLC
  15. * @copyright Copyright (c) 2008, Colorado School of Mines
  16. */
  17. #include "Coupling.h"
  18. #include "FieldPoints.h"
  19. namespace Lemma {
  20. // ==================== FRIEND METHODS =====================
  21. std::ostream &operator << (std::ostream &stream, const Coupling &ob) {
  22. stream << ob.Serialize() << "\n---\n"; // End of doc ---
  23. return stream;
  24. }
  25. // ==================== LIFECYCLE =======================
  26. //--------------------------------------------------------------------------------------
  27. // Class: Coupling
  28. // Method: Coupling
  29. // Description: constructor (locked)
  30. //--------------------------------------------------------------------------------------
  31. Coupling::Coupling (const ctor_key&) : LemmaObject( ) {
  32. } // ----- end of method Coupling::Coupling (constructor) -----
  33. //--------------------------------------------------------------------------------------
  34. // Class: Coupling
  35. // Method: Coupling
  36. // Description: DeSerializing constructor (locked)
  37. //--------------------------------------------------------------------------------------
  38. Coupling::Coupling (const YAML::Node& node, const ctor_key&) : LemmaObject(node) {
  39. } // ----- end of method Coupling::Coupling (constructor) -----
  40. //--------------------------------------------------------------------------------------
  41. // Class: Coupling
  42. // Method: NewSP()
  43. // Description: public constructor returing a shared_ptr
  44. //--------------------------------------------------------------------------------------
  45. std::shared_ptr< Coupling > Coupling::NewSP() {
  46. return std::make_shared< Coupling >( ctor_key() );
  47. }
  48. //--------------------------------------------------------------------------------------
  49. // Class: Coupling
  50. // Method: ~Coupling
  51. // Description: destructor (protected)
  52. //--------------------------------------------------------------------------------------
  53. Coupling::~Coupling () {
  54. } // ----- end of method Coupling::~Coupling (destructor) -----
  55. //--------------------------------------------------------------------------------------
  56. // Class: Coupling
  57. // Method: Serialize
  58. //--------------------------------------------------------------------------------------
  59. YAML::Node Coupling::Serialize ( ) const {
  60. YAML::Node node = LemmaObject::Serialize();
  61. node.SetTag( GetName() );
  62. // Coils Transmitters & Receivers
  63. for ( auto txm : TxRx) {
  64. node[txm.first] = txm.second->Serialize();
  65. }
  66. // LayeredEarthEM
  67. node["SigmaModel"] = SigmaModel->Serialize();
  68. node["tol"] = tol;
  69. node["minLevel"] = minLevel;
  70. node["maxLevel"] = maxLevel;
  71. return node;
  72. } // ----- end of method Coupling::Serialize -----
  73. //--------------------------------------------------------------------------------------
  74. // Class: Coupling
  75. // Method: DeSerialize
  76. //--------------------------------------------------------------------------------------
  77. std::shared_ptr<Coupling> Coupling::DeSerialize ( const YAML::Node& node ) {
  78. if (node.Tag() != "Coupling" ) {
  79. throw DeSerializeTypeMismatch( "Coupling", node.Tag());
  80. }
  81. return std::make_shared< Coupling > ( node, ctor_key() );
  82. } // ----- end of method Coupling::DeSerialize -----
  83. //--------------------------------------------------------------------------------------
  84. // Class: Coupling
  85. // Method: DeSerialize
  86. //--------------------------------------------------------------------------------------
  87. Complex Coupling::Calculate (const std::vector< std::string>& Tx, const std::vector<std::string >& Rx,
  88. bool vtkOutput ) {
  89. static bool first = false; // a little hackish
  90. if (!first) {
  91. // All EM calculations will share same field points
  92. cpoints = FieldPoints::NewSP();
  93. cpoints->SetNumberOfPoints(8);
  94. }
  95. first = true;
  96. for (auto tx : Tx) {
  97. // Set up EMEarth
  98. EMEarths[tx] = EMEarth1D::NewSP();
  99. EMEarths[tx]->AttachWireAntenna(TxRx[tx]);
  100. EMEarths[tx]->AttachLayeredEarthEM(SigmaModel);
  101. EMEarths[tx]->AttachFieldPoints( cpoints );
  102. EMEarths[tx]->SetFieldsToCalculate(H);
  103. // TODO query for method, altough with flat antennae, this is fastest
  104. EMEarths[tx]->SetHankelTransformMethod(ANDERSON801);
  105. EMEarths[tx]->SetTxRxMode(TX);
  106. TxRx[tx]->SetCurrent(1.);
  107. }
  108. for (auto rx : Rx) {
  109. if (EMEarths.count(rx)) {
  110. EMEarths[rx]->SetTxRxMode(TXRX);
  111. } else {
  112. EMEarths[rx] = EMEarth1D::NewSP();
  113. EMEarths[rx]->AttachWireAntenna(TxRx[rx]);
  114. EMEarths[rx]->AttachLayeredEarthEM(SigmaModel);
  115. EMEarths[rx]->AttachFieldPoints( cpoints );
  116. EMEarths[rx]->SetFieldsToCalculate(H);
  117. // TODO query for method, altough with flat antennae, this is fastest
  118. EMEarths[rx]->SetHankelTransformMethod(ANDERSON801);
  119. EMEarths[rx]->SetTxRxMode(RX);
  120. TxRx[rx]->SetCurrent(1.);
  121. }
  122. }
  123. SUM = 0;
  124. IntegrateOnOctreeGrid( vtkOutput );
  125. std::cout << "\nFinished KERNEL\n";
  126. EMEarths.clear();
  127. return SUM;
  128. }
  129. //--------------------------------------------------------------------------------------
  130. // Class: Coupling
  131. // Method: IntegrateOnOctreeGrid
  132. //--------------------------------------------------------------------------------------
  133. void Coupling::IntegrateOnOctreeGrid( bool vtkOutput ) {
  134. static int count = 0;
  135. Vector3r cpos = Origin + Size/2.;
  136. VOLSUM = 0;
  137. nleaves = 0;
  138. if (!vtkOutput) {
  139. EvaluateKids( Size, 0, cpos, Complex(100.));
  140. } else {
  141. #ifdef LEMMAUSEVTK
  142. vtkHyperOctree* oct = vtkHyperOctree::New();
  143. oct->SetDimension(3);
  144. oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
  145. oct->SetSize( Size(0), Size(1), Size(2) );
  146. vtkHyperOctreeCursor* curse = oct->NewCellCursor();
  147. curse->ToRoot();
  148. EvaluateKids2( Size, 0, cpos, Complex(100.0), oct, curse );
  149. // Fill in leaf data
  150. vtkDoubleArray* kr = vtkDoubleArray::New();
  151. kr->SetNumberOfComponents(1);
  152. kr->SetName("Re($K_0$)");
  153. kr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  154. vtkDoubleArray* ki = vtkDoubleArray::New();
  155. ki->SetNumberOfComponents(1);
  156. ki->SetName("Im($K_0$)");
  157. ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  158. vtkDoubleArray* km = vtkDoubleArray::New();
  159. km->SetNumberOfComponents(1);
  160. km->SetName("mod($K_0$)");
  161. km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  162. vtkIntArray* kid = vtkIntArray::New();
  163. kid->SetNumberOfComponents(1);
  164. kid->SetName("ID");
  165. kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  166. vtkIntArray* kerr = vtkIntArray::New();
  167. kerr->SetNumberOfComponents(1);
  168. kerr->SetName("nleaf");
  169. //Real LeafVol(0);
  170. for (auto leaf : LeafDict) {
  171. kr->InsertTuple1( leaf.first, std::real(leaf.second) );
  172. ki->InsertTuple1( leaf.first, std::imag(leaf.second) );
  173. km->InsertTuple1( leaf.first, std::abs(leaf.second) );
  174. kid->InsertTuple1( leaf.first, leaf.first );
  175. //LeafVol += std::real(leaf.second);
  176. }
  177. //std::cout << "\n\nLeafVol=" << LeafVol << std::endl;
  178. for (auto leaf : LeafDictIdx) {
  179. kerr->InsertTuple1( leaf.first, leaf.second );
  180. }
  181. auto kri = oct->GetLeafData()->AddArray(kr);
  182. auto kii = oct->GetLeafData()->AddArray(ki);
  183. auto kmi = oct->GetLeafData()->AddArray(km);
  184. auto kidi = oct->GetLeafData()->AddArray(kid);
  185. auto keri = oct->GetLeafData()->AddArray(kerr);
  186. auto write = vtkXMLHyperOctreeWriter::New();
  187. //write.SetDataModeToAscii()
  188. write->SetInputData(oct);
  189. std::string fname = std::string("octree-couple-") + to_string(count) + std::string(".vto");
  190. write->SetFileName(fname.c_str());
  191. write->Write();
  192. write->Delete();
  193. oct->GetLeafData()->RemoveArray( kri );
  194. oct->GetLeafData()->RemoveArray( kii );
  195. oct->GetLeafData()->RemoveArray( kmi );
  196. oct->GetLeafData()->RemoveArray( kidi );
  197. oct->GetLeafData()->RemoveArray( keri );
  198. kerr->Delete();
  199. kid->Delete();
  200. kr->Delete();
  201. ki->Delete();
  202. km->Delete();
  203. curse->Delete();
  204. oct->Delete();
  205. #else
  206. throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
  207. #endif
  208. }
  209. std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" << Size(0)*Size(1)*Size(2)
  210. << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) << std::endl;
  211. count += 1;
  212. }
  213. //--------------------------------------------------------------------------------------
  214. // Class: Coupling
  215. // Method: f
  216. //--------------------------------------------------------------------------------------
  217. Complex Coupling::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
  218. return volume * ( Ht.dot(Hr) ); // coupling
  219. //return volume * (1.-((Ht+Hr).norm()/(Hr.norm() + Ht.norm()))); // interference
  220. //return volume * std::acos( (Ht.dot(Hr) / (Ht.norm()*Hr.norm())) ); // angle
  221. }
  222. //--------------------------------------------------------------------------------------
  223. // Class: Coupling
  224. // Method: EvaluateKids
  225. //--------------------------------------------------------------------------------------
  226. void Coupling::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
  227. const Complex& parentVal ) {
  228. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  229. std::cout.flush();
  230. // Next level step, interested in one level below
  231. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  232. Vector3r step = size.array() / (Real)(1 << (level+1) );
  233. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  234. Vector3r pos = cpos - step/2.;
  235. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  236. 0, 0, 0,
  237. step[0], 0, 0,
  238. 0, step[1], 0,
  239. step[0], step[1], 0,
  240. 0, 0, step[2],
  241. step[0], 0, step[2],
  242. 0, step[1], step[2],
  243. step[0], step[1], step[2] ).finished();
  244. VectorXcr kvals(8); // individual kernel vals
  245. cpoints->ClearFields();
  246. for (int ichild=0; ichild<8; ++ichild) {
  247. Vector3r cp = pos; // Eigen complains about combining these
  248. cp += posadd.row(ichild);
  249. cpoints->SetLocation( ichild, cp );
  250. }
  251. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  252. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  253. for ( auto EMCalc : EMEarths ) {
  254. EMCalc.second->GetFieldPoints()->ClearFields();
  255. EMCalc.second->CalculateWireAntennaFields();
  256. switch (EMCalc.second->GetTxRxMode()) {
  257. case TX:
  258. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  259. break;
  260. case RX:
  261. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  262. break;
  263. case TXRX:
  264. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  265. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  266. break;
  267. default:
  268. break;
  269. }
  270. }
  271. for (int ichild=0; ichild<8; ++ichild) {
  272. Vector3r cp = pos; // Eigen complains about combining these
  273. cp += posadd.row(ichild);
  274. kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  275. }
  276. Complex ksum = kvals.sum(); // Kernel sum
  277. // Evaluate whether or not furthur splitting is needed
  278. if ( (std::abs(ksum-parentVal) > tol && level < maxLevel) || level < minLevel ) {
  279. // Not a leaf dive further in
  280. for (int ichild=0; ichild<8; ++ichild) {
  281. Vector3r cp = pos; // Eigen complains about combining these
  282. cp += posadd.row(ichild);
  283. EvaluateKids( size, level+1, cp, kvals(ichild) );
  284. }
  285. return; // not leaf
  286. }
  287. // implicit else, is a leaf
  288. SUM += ksum;
  289. VOLSUM += 8.*vol;
  290. nleaves += 1; // could say += 8 just as fairly
  291. return; // is leaf
  292. }
  293. #ifdef LEMMAUSEVTK
  294. //--------------------------------------------------------------------------------------
  295. // Class: Coupling
  296. // Method: EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
  297. //--------------------------------------------------------------------------------------
  298. void Coupling::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
  299. const Complex& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
  300. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  301. std::cout.flush();
  302. // Next level step, interested in one level below
  303. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  304. Vector3r step = size.array() / (Real)(1 << (level+1) );
  305. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  306. Vector3r pos = cpos - step/2.;
  307. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  308. 0, 0, 0,
  309. step[0], 0, 0,
  310. 0, step[1], 0,
  311. step[0], step[1], 0,
  312. 0, 0, step[2],
  313. step[0], 0, step[2],
  314. 0, step[1], step[2],
  315. step[0], step[1], step[2] ).finished();
  316. VectorXcr kvals(8); // individual kernel vals
  317. cpoints->ClearFields();
  318. for (int ichild=0; ichild<8; ++ichild) {
  319. Vector3r cp = pos; // Eigen complains about combining these
  320. cp += posadd.row(ichild);
  321. cpoints->SetLocation( ichild, cp );
  322. }
  323. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  324. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  325. for ( auto EMCalc : EMEarths ) {
  326. //EMCalc->GetFieldPoints()->ClearFields();
  327. EMCalc.second->CalculateWireAntennaFields();
  328. switch (EMCalc.second->GetTxRxMode()) {
  329. case TX:
  330. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  331. break;
  332. case RX:
  333. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  334. break;
  335. case TXRX:
  336. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  337. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  338. break;
  339. default:
  340. break;
  341. }
  342. }
  343. for (int ichild=0; ichild<8; ++ichild) {
  344. Vector3r cp = pos; // Eigen complains about combining these
  345. cp += posadd.row(ichild);
  346. kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  347. }
  348. Complex ksum = kvals.sum(); // Kernel sum
  349. // Evaluate whether or not furthur splitting is needed
  350. if ( (std::abs(ksum-parentVal) > tol && level < maxLevel) || level < minLevel ) {
  351. oct->SubdivideLeaf(curse);
  352. for (int ichild=0; ichild<8; ++ichild) {
  353. curse->ToChild(ichild);
  354. Vector3r cp = pos; // Eigen complains about combining these
  355. cp += posadd.row(ichild);
  356. /* Test for position via alternative means */
  357. /*
  358. Real p[3];
  359. GetPosition(curse, p);
  360. if ( (Vector3r(p) - cp).norm() > 1e-8 ) {
  361. std::cout << "ERROR @ nleaves" << nleaves << "\n" << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
  362. << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
  363. throw std::runtime_error("doom");
  364. }
  365. */
  366. /* End of position test */
  367. EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
  368. curse->ToParent();
  369. }
  370. return; // not a leaf
  371. }
  372. LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
  373. LeafDictIdx[curse->GetLeafId()] = nleaves;
  374. SUM += ksum;
  375. VOLSUM += 8*vol;
  376. nleaves += 1;
  377. return; // is a leaf
  378. }
  379. //--------------------------------------------------------------------------------------
  380. // Class: Coupling
  381. // Method: GetPosition
  382. //--------------------------------------------------------------------------------------
  383. void Coupling::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
  384. Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
  385. //step = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
  386. p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
  387. p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
  388. p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
  389. }
  390. #endif
  391. } // ---- end of namespace Lemma ----
  392. /* vim: set tabstop=4 expandtab */
  393. /* vim: set filetype=cpp */