Surface NMR forward modelling
Du kan inte välja fler än 25 ämnen Ämnen måste starta med en bokstav eller siffra, kan innehålla bindestreck ('-') och vara max 35 tecken långa.

KernelV0.cpp 31KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693
  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 "MerlinConfig.h"
  18. #include "KernelV0.h"
  19. #include "FieldPoints.h"
  20. namespace Lemma {
  21. // ==================== FRIEND METHODS =====================
  22. std::ostream &operator << (std::ostream &stream, const KernelV0 &ob) {
  23. stream << ob.Serialize() << "\n"; // End of doc ---
  24. return stream;
  25. }
  26. // ==================== LIFECYCLE =======================
  27. //--------------------------------------------------------------------------------------
  28. // Class: KernelV0
  29. // Method: KernelV0
  30. // Description: constructor (locked)
  31. //--------------------------------------------------------------------------------------
  32. KernelV0::KernelV0 (const ctor_key& key) : MerlinObject( key ) {
  33. } // ----- end of method KernelV0::KernelV0 (constructor) -----
  34. //--------------------------------------------------------------------------------------
  35. // Class: KernelV0
  36. // Method: KernelV0
  37. // Description: DeSerializing constructor (locked)
  38. //--------------------------------------------------------------------------------------
  39. KernelV0::KernelV0 (const YAML::Node& node, const ctor_key& key) : MerlinObject(node, key) {
  40. //node["PulseType"] = "FID";
  41. Larmor = node["Larmor"].as<Real>();
  42. Temperature = node["Temperature"].as<Real>();
  43. tol = node["tol"].as<Real>();
  44. minLevel = node["minLevel"].as<int>();
  45. maxLevel = node["maxLevel"].as<int>();
  46. Interfaces = node["Interfaces"].as<VectorXr>();
  47. Size = node["IntegrationSize"].as<Vector3r>();
  48. Origin = node["IntegrationOrigin"].as<Vector3r>();
  49. if (node["AlignWithAkvoData"]) {
  50. // Match pulse info with dataset
  51. AlignWithAkvoDataset( YAML::LoadFile( node["AlignWithAkvoData"].as<std::string>()));
  52. } else {
  53. // Read Pulse info direct from Kernel file
  54. PulseI = node["PulseI"].as<VectorXr>();
  55. Taup = node["Taup"].as<Real>();
  56. }
  57. if (node["SigmaModel"]) {
  58. if (node["SigmaModel"].Tag() == "LayeredEarthEM") {
  59. SigmaModel = LayeredEarthEM::DeSerialize(node["SigmaModel"]);
  60. } else {
  61. SigmaModel = LayeredEarthEM::DeSerialize( YAML::LoadFile( node["SigmaModel"].as<std::string>() ));
  62. }
  63. }
  64. if (node["Coils"]) {
  65. for ( auto coil : node["Coils"] ) {
  66. if ( coil.second.Tag() == "PolygonalWireAntenna" ) {
  67. TxRx[ coil.first.as<std::string>() ] = PolygonalWireAntenna::DeSerialize( coil.second );
  68. } else {
  69. TxRx[ coil.first.as<std::string>() ] =
  70. PolygonalWireAntenna::DeSerialize( YAML::LoadFile(coil.second.as<std::string>()) );
  71. }
  72. }
  73. }
  74. if (node["K0"]) {
  75. Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size() ).array() + 1.;
  76. for ( int ilay=0; ilay<Interfaces.size()-1; ++ilay ) {
  77. Kern.row(ilay) = node["K0"]["layer-" + to_string(ilay) ].as<VectorXcr>();
  78. }
  79. }
  80. } // ----- end of method KernelV0::KernelV0 (constructor) -----
  81. //--------------------------------------------------------------------------------------
  82. // Class: KernelV0
  83. // Method: NewSP()
  84. // Description: public constructor returing a shared_ptr
  85. //--------------------------------------------------------------------------------------
  86. std::shared_ptr< KernelV0 > KernelV0::NewSP() {
  87. return std::make_shared< KernelV0 >( ctor_key() );
  88. }
  89. //--------------------------------------------------------------------------------------
  90. // Class: KernelV0
  91. // Method: ~KernelV0
  92. // Description: destructor (protected)
  93. //--------------------------------------------------------------------------------------
  94. KernelV0::~KernelV0 () {
  95. } // ----- end of method KernelV0::~KernelV0 (destructor) -----
  96. //--------------------------------------------------------------------------------------
  97. // Class: KernelV0
  98. // Method: Serialize
  99. //--------------------------------------------------------------------------------------
  100. YAML::Node KernelV0::Serialize ( ) const {
  101. YAML::Node node = MerlinObject::Serialize();
  102. node.SetTag( GetName() );
  103. // Coils Transmitters & Receivers
  104. if (!TxRx.empty()) {
  105. for ( auto txm : TxRx) {
  106. node["Coils"][txm.first] = txm.second->Serialize();
  107. }
  108. }
  109. // LayeredEarthEM
  110. if (SigmaModel != nullptr) {
  111. node["SigmaModel"] = SigmaModel->Serialize();
  112. }
  113. node["PulseType"] = "FID";
  114. node["Larmor"] = Larmor;
  115. node["Temperature"] = Temperature;
  116. node["tol"] = tol;
  117. node["minLevel"] = minLevel;
  118. node["maxLevel"] = maxLevel;
  119. node["Taup"] = Taup;
  120. node["PulseI"] = PulseI;
  121. node["Interfaces"] = Interfaces;
  122. node["IntegrationSize"] = Size;
  123. node["IntegrationOrigin"] = Origin;
  124. // TODO, use better matrix encapulation
  125. if (Kern.array().abs().any() > 1e-16) {
  126. for ( int ilay=0; ilay<Interfaces.size()-1; ++ilay ) {
  127. node["K0"]["layer-" + to_string(ilay) ] = static_cast<VectorXcr>(Kern.row(ilay));
  128. }
  129. }
  130. return node;
  131. } // ----- end of method KernelV0::Serialize -----
  132. //--------------------------------------------------------------------------------------
  133. // Class: KernelV0
  134. // Method: DeSerialize
  135. //--------------------------------------------------------------------------------------
  136. std::shared_ptr<KernelV0> KernelV0::DeSerialize ( const YAML::Node& node ) {
  137. if (node.Tag() != "KernelV0" ) {
  138. throw DeSerializeTypeMismatch( "KernelV0", node.Tag());
  139. }
  140. return std::make_shared< KernelV0 > ( node, ctor_key() );
  141. } // ----- end of method KernelV0::DeSerialize -----
  142. //--------------------------------------------------------------------------------------
  143. // Class: KernelV0
  144. // Method: AlignWithAkvoDataset
  145. //--------------------------------------------------------------------------------------
  146. void KernelV0::AlignWithAkvoDataset( const YAML::Node& node ) {
  147. if (node["processed"].as<std::string>().substr(0,4) == "Akvo") {
  148. std::cout << "Akvo file read\n";
  149. std::cout << node["processed"] << std::endl;
  150. }
  151. if (node["pulseType"].as<std::string>() == "FID") {
  152. std::cout << "FID pulse detected" << std::endl;
  153. PulseI = node["Pulses"]["Pulse 1"]["current"].as<VectorXr>();
  154. this->SetPulseDuration( node["pulseLength"][0].as<double>() );
  155. } else {
  156. std::cerr << "Pulse Type " << node["PulseType"] << "is not supported\n";
  157. }
  158. std::cout << "Finished with Akvo file read" << std::endl;
  159. }
  160. //--------------------------------------------------------------------------------------
  161. // Class: KernelV0
  162. // Method: DeSerialize
  163. //--------------------------------------------------------------------------------------
  164. void KernelV0::CalculateK0 (const std::vector< std::string>& Tx,
  165. const std::vector<std::string >& Rx, bool vtkOutput ) {
  166. // Set up
  167. Larmor = SigmaModel->GetMagneticFieldMagnitude()*GAMMA; // in rad 2246.*2.*PI;
  168. // All EM calculations will share same field points
  169. cpoints = FieldPoints::NewSP();
  170. cpoints->SetNumberOfPoints(8);
  171. for (auto tx : Tx) {
  172. // Set up EMEarth
  173. EMEarths[tx] = EMEarth1D::NewSP();
  174. EMEarths[tx]->AttachWireAntenna(TxRx[tx]);
  175. EMEarths[tx]->AttachLayeredEarthEM(SigmaModel);
  176. EMEarths[tx]->AttachFieldPoints( cpoints );
  177. EMEarths[tx]->SetFieldsToCalculate(H);
  178. // TODO query for method, altough with flat antennae, this is fastest
  179. //EMEarths[tx]->SetHankelTransformMethod(FHTKEY201);
  180. EMEarths[tx]->SetHankelTransformMethod(ANDERSON801);
  181. EMEarths[tx]->SetTxRxMode(TX);
  182. TxRx[tx]->SetCurrent(1.);
  183. }
  184. for (auto rx : Rx) {
  185. if (EMEarths.count(rx)) {
  186. EMEarths[rx]->SetTxRxMode(TXRX);
  187. } else {
  188. EMEarths[rx] = EMEarth1D::NewSP();
  189. EMEarths[rx]->AttachWireAntenna(TxRx[rx]);
  190. EMEarths[rx]->AttachLayeredEarthEM(SigmaModel);
  191. EMEarths[rx]->AttachFieldPoints( cpoints );
  192. EMEarths[rx]->SetFieldsToCalculate(H);
  193. // TODO query for method, altough with flat antennae, this is fastest
  194. //EMEarths[rx]->SetHankelTransformMethod(FHTKEY201);
  195. EMEarths[rx]->SetHankelTransformMethod(ANDERSON801);
  196. EMEarths[rx]->SetTxRxMode(RX);
  197. TxRx[rx]->SetCurrent(1.);
  198. }
  199. }
  200. std::cout << "Calculating K0 kernel\n";
  201. Kern = MatrixXcr::Zero( Interfaces.size()-1, PulseI.size() );
  202. for (ilay=0; ilay<Interfaces.size()-1; ++ilay) {
  203. std::cout << "Layer " << ilay << "\tfrom " << Interfaces(ilay) <<" to "
  204. << Interfaces(ilay+1) << std::endl;
  205. Size(2) = Interfaces(ilay+1) - Interfaces(ilay);
  206. Origin(2) = Interfaces(ilay);
  207. IntegrateOnOctreeGrid( vtkOutput );
  208. }
  209. std::cout << "\nFinished KERNEL\n";
  210. }
  211. //--------------------------------------------------------------------------------------
  212. // Class: KernelV0
  213. // Method: IntegrateOnOctreeGrid
  214. //--------------------------------------------------------------------------------------
  215. void KernelV0::IntegrateOnOctreeGrid( bool vtkOutput) {
  216. Vector3r cpos = Origin + Size/2.;
  217. VOLSUM = 0;
  218. nleaves = 0;
  219. if (!vtkOutput) {
  220. EvaluateKids( Size, 0, cpos, VectorXcr::Ones(PulseI.size()) );
  221. } else {
  222. #ifdef LEMMAUSEVTK
  223. vtkHyperOctree* oct = vtkHyperOctree::New();
  224. oct->SetDimension(3);
  225. oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
  226. oct->SetSize( Size(0), Size(1), Size(2) );
  227. vtkHyperOctreeCursor* curse = oct->NewCellCursor();
  228. curse->ToRoot();
  229. EvaluateKids2( Size, 0, cpos, VectorXcr::Ones(PulseI.size()), oct, curse );
  230. for (int iq=0; iq<PulseI.size(); ++iq) {
  231. // Fill in leaf data
  232. vtkDoubleArray* kr = vtkDoubleArray::New();
  233. kr->SetNumberOfComponents(1);
  234. kr->SetName("Re($\\mathcal{K}_0$)");
  235. kr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  236. vtkDoubleArray* ki = vtkDoubleArray::New();
  237. ki->SetNumberOfComponents(1);
  238. ki->SetName("Im($\\mathcal{K}_0$)");
  239. ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  240. vtkDoubleArray* km = vtkDoubleArray::New();
  241. km->SetNumberOfComponents(1);
  242. km->SetName("mod($\\mathcal{K}_0$)");
  243. km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  244. vtkIntArray* kid = vtkIntArray::New();
  245. kid->SetNumberOfComponents(1);
  246. kid->SetName("ID");
  247. kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  248. vtkIntArray* kerr = vtkIntArray::New();
  249. kerr->SetNumberOfComponents(1);
  250. kerr->SetName("err");
  251. kerr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  252. // Ht field
  253. vtkDoubleArray* htr = vtkDoubleArray::New();
  254. htr->SetNumberOfComponents(3);
  255. htr->SetName("Re($\\mathbf{\\mathcal{H}}_T$)");
  256. htr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  257. vtkDoubleArray* hti = vtkDoubleArray::New();
  258. hti->SetNumberOfComponents(3);
  259. hti->SetName("Im($\\mathbf{\\mathcal{H}}_T$)");
  260. hti->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  261. // Hr field
  262. vtkDoubleArray* hrr = vtkDoubleArray::New();
  263. hrr->SetNumberOfComponents(3);
  264. hrr->SetName("Re($\\mathbf{\\mathcal{H}}_R$)");
  265. hrr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  266. vtkDoubleArray* hri = vtkDoubleArray::New();
  267. hri->SetNumberOfComponents(3);
  268. hri->SetName("Im($\\mathbf{\\mathcal{H}}_R$)");
  269. hri->SetNumberOfTuples( oct->GetNumberOfLeaves() );
  270. //Real LeafVol(0);
  271. for (auto leaf : LeafDict) {
  272. kr->InsertTuple1( leaf.first, std::real(leaf.second(iq)) );
  273. ki->InsertTuple1( leaf.first, std::imag(leaf.second(iq)) );
  274. km->InsertTuple1( leaf.first, std::abs(leaf.second(iq)) );
  275. kid->InsertTuple1( leaf.first, leaf.first );
  276. //LeafVol += std::real(leaf.second);
  277. }
  278. for (auto leaf : LeafHt ) {
  279. htr->InsertTuple( leaf.first, leaf.second.real().data() );
  280. hti->InsertTuple( leaf.first, leaf.second.imag().data() );
  281. }
  282. for (auto leaf : LeafHr ) {
  283. hrr->InsertTuple( leaf.first, leaf.second.real().data() );
  284. hri->InsertTuple( leaf.first, leaf.second.imag().data() );
  285. }
  286. for (auto leaf : LeafDictIdx) {
  287. kerr->InsertTuple1( leaf.first, leaf.second );
  288. }
  289. auto kri = oct->GetLeafData()->AddArray(kr);
  290. auto kii = oct->GetLeafData()->AddArray(ki);
  291. auto kmi = oct->GetLeafData()->AddArray(km);
  292. auto kidi = oct->GetLeafData()->AddArray(kid);
  293. auto keri = oct->GetLeafData()->AddArray(kerr);
  294. auto khtr = oct->GetLeafData()->AddArray(htr);
  295. auto khti = oct->GetLeafData()->AddArray(hti);
  296. auto khrr = oct->GetLeafData()->AddArray(hrr);
  297. auto khri = oct->GetLeafData()->AddArray(hri);
  298. auto write = vtkXMLHyperOctreeWriter::New();
  299. //write.SetDataModeToAscii()
  300. write->SetInputData(oct);
  301. std::string fname = std::string("octree-") + to_string(ilay)
  302. + std::string("-") + to_string(iq) + std::string(".vto");
  303. write->SetFileName(fname.c_str());
  304. write->Write();
  305. write->Delete();
  306. oct->GetLeafData()->RemoveArray( kri );
  307. oct->GetLeafData()->RemoveArray( kii );
  308. oct->GetLeafData()->RemoveArray( kmi );
  309. oct->GetLeafData()->RemoveArray( kidi );
  310. oct->GetLeafData()->RemoveArray( keri );
  311. oct->GetLeafData()->RemoveArray( khtr );
  312. oct->GetLeafData()->RemoveArray( khti );
  313. oct->GetLeafData()->RemoveArray( khrr );
  314. oct->GetLeafData()->RemoveArray( khri );
  315. kerr->Delete();
  316. kid->Delete();
  317. kr->Delete();
  318. ki->Delete();
  319. km->Delete();
  320. htr->Delete();
  321. hti->Delete();
  322. hrr->Delete();
  323. hri->Delete();
  324. }
  325. curse->Delete();
  326. oct->Delete();
  327. #else
  328. throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
  329. #endif
  330. }
  331. std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" << Size(0)*Size(1)*Size(2)
  332. << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) << std::endl;
  333. }
  334. //--------------------------------------------------------------------------------------
  335. // Class: KernelV0
  336. // Method: f
  337. //--------------------------------------------------------------------------------------
  338. VectorXcr KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
  339. // Compute the elliptic fields
  340. Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
  341. Vector3r B0 = SigmaModel->GetMagneticField();
  342. // Elliptic representation
  343. EllipticB EBT = EllipticFieldRep(MU0*Ht, B0hat);
  344. EllipticB EBR = EllipticFieldRep(MU0*Hr, B0hat);
  345. // Compute Mn0
  346. Vector3r Mn0 = ComputeMn0(1.0, B0);
  347. Real Mn0Abs = Mn0.norm();
  348. //std::cout << "Mn0\t" << Mn0.transpose() << std::endl;
  349. // Compute phase delay
  350. // TODO add transmiiter current phase and delay induced apparent time phase!
  351. Complex PhaseTerm = EBR.bhat.dot(EBT.bhat) + Complex(0, (B0hat.dot(EBR.bhat.cross(EBT.bhat))));
  352. Complex ejztr = std::exp(Complex(0, EBR.zeta + EBT.zeta));
  353. // Calcuate vector of all responses
  354. VectorXcr F = VectorXcr::Zero( PulseI.size() );
  355. for (int iq=0; iq<PulseI.size(); ++iq) {
  356. // Compute the tipping angle
  357. Real sintheta = std::sin(0.5*GAMMA*PulseI(iq)*Taup*(EBT.alpha-EBT.beta));
  358. F(iq) = -volume*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
  359. }
  360. return F;
  361. }
  362. // //--------------------------------------------------------------------------------------
  363. // // Class: KernelV0
  364. // // Method: ComputeV0Cell
  365. // //--------------------------------------------------------------------------------------
  366. // Complex KernelV0::ComputeV0Cell(const EllipticB& EBT, const EllipticB& EBR,
  367. // const Real& sintheta, const Real& phase, const Real& Mn0Abs,
  368. // const Real& vol) {
  369. // // earth response of receiver adjoint field
  370. // Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
  371. // Complex ejztr = std::exp(Complex(0, EBR.zeta + EBT.zeta));
  372. // Complex PhaseTerm = EBR.bhat.dot(EBT.bhat) + (B0hat.dot(EBR.bhat.cross(EBT.bhat) ));
  373. // return -vol*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
  374. // }
  375. //--------------------------------------------------------------------------------------
  376. // Class: KernelV0
  377. // Method: ComputeV0Cell
  378. //--------------------------------------------------------------------------------------
  379. Vector3r KernelV0::ComputeMn0(const Real& Porosity, const Vector3r& B0) {
  380. Real chi_n = NH2O*((GAMMA*GAMMA*HBAR*HBAR)/(4.*KB*Temperature));
  381. return chi_n*Porosity*B0;
  382. }
  383. //--------------------------------------------------------------------------------------
  384. // Class: KernelV0
  385. // Method: ComputeV0Cell
  386. //--------------------------------------------------------------------------------------
  387. EllipticB KernelV0::EllipticFieldRep (const Vector3cr& B, const Vector3r& B0hat) {
  388. // This all follows Weichman et al., 2000.
  389. // There are some numerical stability issues that arise when the two terms in the beta
  390. // formulation are nearly equivalent. The current formulation will result in a null-valued
  391. // beta, or can underflow. However, this does not entirely recreate the true value of B perp.
  392. // Error is checked to be below 1%, but reformulating for numeric stability may be welcome
  393. EllipticB ElipB = EllipticB();
  394. Vector3cr Bperp = B - B0hat.dot(B)*B0hat;
  395. Real BperpNorm = Bperp.norm();
  396. // These two are equivalent
  397. //Complex Bp2 = Bperp.transpose() * Bperp;
  398. Complex Bp2 = Bperp.conjugate().dot(Bperp);
  399. VectorXcr iB0 = Complex(0,1)*B0hat.cast<Complex>().array();
  400. ElipB.eizt = std::sqrt(Bp2 / std::abs(Bp2));
  401. ElipB.alpha = INVSQRT2*std::sqrt(BperpNorm*BperpNorm + std::abs(Bp2));
  402. //ElipB.beta = std::copysign(1, std::real(iB0.dot( Bperp.cross(Bperp.conjugate())) )) *
  403. ElipB.beta = sgn( std::real(iB0.dot( Bperp.cross(Bperp.conjugate())) )) *
  404. (INVSQRT2*std::sqrt(BperpNorm*BperpNorm - std::abs(Bp2)));
  405. // Correct underflow in beta calculation
  406. // could use cerrno instead...
  407. // http://en.cppreference.com/w/cpp/numeric/math/sqrt
  408. if (ElipB.beta != ElipB.beta) ElipB.beta = 0;
  409. ElipB.bhat = ((Real)1./ElipB.alpha)*(((Real)1./ElipB.eizt)*Bperp.array()).real().array();
  410. ElipB.bhatp = B0hat.cross(ElipB.bhat);
  411. ElipB.zeta = std::real(std::log(ElipB.eizt)/Complex(0,1));
  412. /* as an error check decomposed field - computed actual */
  413. // Vector3cr Bperp2 = ElipB.eizt * (ElipB.alpha * ElipB.bhat
  414. // + (Complex(0,1) * ElipB.beta * ElipB.bhatp) );
  415. // ElipB.err = (Bperp-Bperp2).norm();
  416. // if (ElipB.err > .01*Bperp.norm() ) { // 1% error
  417. // std::cout << "Elip error\n";
  418. // Real Beta2 = sgn( std::real(iB0.dot( Bperp.cross(Bperp.conjugate())) )) *
  419. // (INVSQRT2*std::sqrt(BperpNorm*BperpNorm - std::abs(Bp2)));
  420. // Vector3cr Bperp3 = ElipB.eizt * (ElipB.alpha * ElipB.bhat
  421. // + (Complex(0,1) * Beta2 * ElipB.bhatp) );
  422. // std::cout << "Beta term0\t" << (INVSQRT2*std::sqrt(BperpNorm*BperpNorm - std::abs(Bp2))) << std::endl;
  423. // std::cout << "Beta term1\t" << BperpNorm*BperpNorm << "\t" << std::abs(Bp2) << std::endl;
  424. // std::cout << "Beta \t" << ElipB.beta << std::endl;
  425. // std::cout << "Beta2 \t" << Beta2 << std::endl;
  426. // std::cout << "Bperp \t" << Bperp.transpose() << std::endl;
  427. // std::cout << "Bperp2\t" << Bperp2.transpose() << std::endl;
  428. // std::cout << "Bperp3\t" << Bperp3.transpose() << std::endl;
  429. // std::cout << "err \t" << ElipB.err << std::endl;
  430. // }
  431. return ElipB;
  432. }
  433. //--------------------------------------------------------------------------------------
  434. // Class: KernelV0
  435. // Method: EvaluateKids
  436. //--------------------------------------------------------------------------------------
  437. void KernelV0::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
  438. const VectorXcr& parentVal ) {
  439. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  440. //std::cout.flush();
  441. // Next level step, interested in one level below
  442. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  443. Vector3r step = size.array() / (Real)(1 << (level+1) );
  444. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  445. Vector3r pos = cpos - step/2.;
  446. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  447. 0, 0, 0,
  448. step[0], 0, 0,
  449. 0, step[1], 0,
  450. step[0], step[1], 0,
  451. 0, 0, step[2],
  452. step[0], 0, step[2],
  453. 0, step[1], step[2],
  454. step[0], step[1], step[2] ).finished();
  455. cpoints->ClearFields();
  456. for (int ichild=0; ichild<8; ++ichild) {
  457. Vector3r cp = pos; // Eigen complains about combining these
  458. cp += posadd.row(ichild);
  459. cpoints->SetLocation( ichild, cp );
  460. }
  461. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  462. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  463. for ( auto EMCalc : EMEarths ) {
  464. EMCalc.second->GetFieldPoints()->ClearFields();
  465. EMCalc.second->CalculateWireAntennaFields();
  466. switch (EMCalc.second->GetTxRxMode()) {
  467. case TX:
  468. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  469. break;
  470. case RX:
  471. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  472. break;
  473. case TXRX:
  474. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  475. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  476. break;
  477. default:
  478. break;
  479. }
  480. }
  481. MatrixXcr kvals(8, PulseI.size()); // individual kernel vals
  482. for (int ichild=0; ichild<8; ++ichild) {
  483. Vector3r cp = pos; // Eigen complains about combining these
  484. cp += posadd.row(ichild);
  485. kvals.row(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  486. }
  487. VectorXcr ksum = kvals.colwise().sum(); // Kernel sum
  488. // Evaluate whether or not furthur splitting is needed
  489. if ( (((ksum - parentVal).array().abs() > tol).any() && level<maxLevel) || level < minLevel ) {
  490. // Not a leaf dive further in
  491. for (int ichild=0; ichild<8; ++ichild) {
  492. Vector3r cp = pos; // Eigen complains about combining these
  493. cp += posadd.row(ichild);
  494. EvaluateKids( size, level+1, cp, kvals.row(ichild) );
  495. }
  496. return; // not leaf
  497. }
  498. // implicit else, is a leaf
  499. Kern.row(ilay) += ksum;
  500. VOLSUM += 8.*vol;
  501. nleaves += 8; // reflects the number of kernel evaluations
  502. return; // is leaf
  503. }
  504. #ifdef LEMMAUSEVTK
  505. //--------------------------------------------------------------------------------------
  506. // Class: KernelV0
  507. // Method: EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
  508. //--------------------------------------------------------------------------------------
  509. void KernelV0::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
  510. const VectorXcr& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
  511. std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
  512. //std::cout.flush();
  513. // Next level step, interested in one level below
  514. // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
  515. Vector3r step = size.array() / (Real)(1 << (level+1) );
  516. Real vol = (step(0)*step(1)*step(2)); // volume of each child
  517. Vector3r pos = cpos - step/2.;
  518. Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
  519. 0, 0, 0,
  520. step[0], 0, 0,
  521. 0, step[1], 0,
  522. step[0], step[1], 0,
  523. 0, 0, step[2],
  524. step[0], 0, step[2],
  525. 0, step[1], step[2],
  526. step[0], step[1], step[2] ).finished();
  527. MatrixXcr kvals(8, PulseI.size()); // individual kernel vals
  528. cpoints->ClearFields();
  529. for (int ichild=0; ichild<8; ++ichild) {
  530. Vector3r cp = pos; // Eigen complains about combining these
  531. cp += posadd.row(ichild);
  532. cpoints->SetLocation( ichild, cp );
  533. }
  534. Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
  535. Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
  536. for ( auto EMCalc : EMEarths ) {
  537. //EMCalc->GetFieldPoints()->ClearFields();
  538. EMCalc.second->CalculateWireAntennaFields();
  539. switch (EMCalc.second->GetTxRxMode()) {
  540. case TX:
  541. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  542. break;
  543. case RX:
  544. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  545. break;
  546. case TXRX:
  547. Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
  548. Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
  549. break;
  550. default:
  551. break;
  552. }
  553. }
  554. for (int ichild=0; ichild<8; ++ichild) {
  555. Vector3r cp = pos; // Eigen complains about combining these
  556. cp += posadd.row(ichild);
  557. kvals.row(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
  558. }
  559. VectorXcr ksum = kvals.colwise().sum(); // Kernel sum
  560. // Evaluate whether or not furthur splitting is needed
  561. if ( (((ksum - parentVal).array().abs() > tol).any() && level<maxLevel) || level < minLevel ) {
  562. oct->SubdivideLeaf(curse);
  563. for (int ichild=0; ichild<8; ++ichild) {
  564. curse->ToChild(ichild);
  565. Vector3r cp = pos; // Eigen complains about combining these
  566. cp += posadd.row(ichild);
  567. /* Test for position via alternative means */
  568. /*
  569. Real p[3];
  570. GetPosition(curse, p);
  571. if ( (Vector3r(p) - cp).norm() > 1e-8 ) {
  572. std::cout << "ERROR @ nleaves" << nleaves << "\n" << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
  573. << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
  574. throw std::runtime_error("doom");
  575. }
  576. */
  577. /* End of position test */
  578. EvaluateKids2( size, level+1, cp, kvals.row(ichild), oct, curse );
  579. curse->ToParent();
  580. }
  581. return; // not a leaf
  582. }
  583. /* just stuff with sum of the kids and don't subdivide */
  584. /*
  585. LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
  586. LeafDictIdx[curse->GetLeafId()] = nleaves;
  587. */
  588. /* Alternatively, subdivide the VTK octree here and stuff the children to make better
  589. * visuals, but also 8x the storage...
  590. */
  591. oct->SubdivideLeaf(curse);
  592. for (int ichild=0; ichild<8; ++ichild) {
  593. curse->ToChild(ichild);
  594. LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
  595. LeafHt[curse->GetLeafId()] = Ht.col(ichild);
  596. LeafHr[curse->GetLeafId()] = Hr.col(ichild);
  597. LeafDictIdx[curse->GetLeafId()] = nleaves;
  598. curse->ToParent();
  599. }
  600. Kern.row(ilay) += ksum;
  601. VOLSUM += 8*vol;
  602. nleaves += 8; // good reason to say 1 or 8 here...8 sounds better and reflects kernel evaluations
  603. return; // is a leaf
  604. }
  605. //--------------------------------------------------------------------------------------
  606. // Class: KernelV0
  607. // Method: GetPosition
  608. //--------------------------------------------------------------------------------------
  609. void KernelV0::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
  610. Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
  611. //step = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
  612. p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
  613. p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
  614. p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
  615. }
  616. #endif
  617. } // ---- end of namespace Lemma ----
  618. /* vim: set tabstop=4 expandtab */
  619. /* vim: set filetype=cpp */