Surface NMR forward modelling
Du kannst nicht mehr als 25 Themen auswählen Themen müssen mit entweder einem Buchstaben oder einer Ziffer beginnen. Sie können Bindestriche („-“) enthalten und bis zu 35 Zeichen lang sein.

KernelV0.h 15KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416
  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:34 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. #ifndef KERNELV0_INC
  18. #define KERNELV0_INC
  19. #pragma once
  20. #include "MerlinObject.h"
  21. #include "LayeredEarthEM.h"
  22. #include "PolygonalWireAntenna.h"
  23. #include "EMEarth1D.h"
  24. #ifdef LEMMAUSEVTK
  25. //#include "vtkHyperOctree.h"
  26. //#include "vtkHyperOctreeCursor.h"
  27. //#include "vtkXMLHyperOctreeWriter.h"
  28. //#include "vtkXMLHyperOctreeWriter.h"
  29. #include "vtkCellData.h"
  30. #include "vtkPointData.h"
  31. #include "vtkHyperTree.h"
  32. #include "vtkHyperTree.h"
  33. #include "vtkHyperTreeGrid.h"
  34. #include "vtkXMLHyperTreeGridWriter.h"
  35. #include "vtkHyperTreeCursor.h" // not in VTK 8.90
  36. //#include "vtkHyperTreeGridLevelEntry.h" VTK 9
  37. #include "vtkDoubleArray.h"
  38. #endif
  39. namespace Lemma {
  40. // Holds the elliptic field construction of Bperp
  41. // commented out variables are for error checking
  42. struct EllipticB {
  43. Real alpha;
  44. Real beta;
  45. Real zeta;
  46. // Real err;
  47. Complex eizt;
  48. // Complex BperpdotB;
  49. Vector3r bhat;
  50. Vector3r bhatp;
  51. // Vector3cr Bperp;
  52. };
  53. template <typename T> int sgn(T val) {
  54. return (val > T(0)) - (val < T(0));
  55. }
  56. /**
  57. * \ingroup Merlin
  58. * \brief Calculated the initial amplitude imaging kernel of a sNMR experiment
  59. * \details This class calculates the imaging kernel for a free induction decay
  60. * pulse. The methodology follows from Weichman et al., 2000.
  61. */
  62. class KernelV0 : public MerlinObject {
  63. friend std::ostream &operator<<(std::ostream &stream, const KernelV0 &ob);
  64. protected:
  65. /*
  66. * This key is used to lock the constructor. It is protected so that inhereted
  67. * classes also have the key to contruct their base class.
  68. */
  69. public:
  70. // ==================== LIFECYCLE =======================
  71. /**
  72. * Default constructor.
  73. * @note This method is locked, and cannot be called directly.
  74. * The reason that the method is public is to enable the use
  75. * of make_shared whilst enforcing the use of shared_ptr,
  76. * in c++-17, this curiosity may be resolved.
  77. * @see KernelV0::NewSP
  78. */
  79. explicit KernelV0 ( const ctor_key& );
  80. /**
  81. * DeSerializing constructor.
  82. * @note This method is locked, and cannot be called directly.
  83. * The reason that the method is public is to enable the use
  84. * of make_shared whilst enforcing the use of shared_ptr,
  85. * in c++-17, this curiosity may be resolved.
  86. * @see KernelV0::DeSerialize
  87. */
  88. KernelV0 ( const YAML::Node& node, const ctor_key& );
  89. /**
  90. * Default destructor.
  91. * @note This method should never be called due to the mandated
  92. * use of smart pointers. It is necessary to keep the method
  93. * public in order to allow for the use of the more efficient
  94. * make_shared constructor.
  95. */
  96. virtual ~KernelV0 ();
  97. /**
  98. * Uses YAML to serialize this object.
  99. * @return a YAML::Node
  100. * @see KernelV0::DeSerialize
  101. */
  102. virtual YAML::Node Serialize() const;
  103. /**
  104. * Factory method for generating concrete class.
  105. * @return a std::shared_ptr of type KernelV0
  106. */
  107. static std::shared_ptr< KernelV0 > NewSP();
  108. /**
  109. * Constructs an KernelV0 object from a YAML::Node.
  110. * @see KernelV0::Serialize
  111. */
  112. static std::shared_ptr<KernelV0> DeSerialize(const YAML::Node& node);
  113. /**
  114. * Constructs an object from a string representation of a YAML::Node. This is primarily
  115. * used in Python wrapping
  116. */
  117. static std::shared_ptr<KernelV0> DeSerialize( const std::string& node ) {
  118. return KernelV0::DeSerialize(YAML::LoadFile(node));
  119. }
  120. // ==================== OPERATORS =======================
  121. // ==================== OPERATIONS =======================
  122. /**
  123. * Calculates a single imaging kernel, however, phased arrays are supported
  124. * so that more than one transmitter and/or receiver can be specified.
  125. * @param[in] tx is the list of transmitters to use for a kernel, use the same labels as
  126. * used in PushCoil.
  127. * @param[in] rx is the list of receivers to use for a kernel, use the same labels as
  128. * used in PushCoil. @see PushCoil
  129. * @param[in] vtkOutput generates a VTK hyperoctree file as well, useful for visualization.
  130. * requires compilation of Lemma with VTK. The VTK files can become very large.
  131. */
  132. void CalculateK0 (const std::vector< std::string >& tx, const std::vector< std::string >& rx,
  133. bool vtkOutput );
  134. //bool vtkOutput=false );
  135. /**
  136. * Aligns the kernel pulse settings with an Akvo Processed dataset.
  137. */
  138. void AlignWithAkvoDataset( const YAML::Node& node ) ;
  139. /**
  140. * Aligns the kernel pulse settings with an Akvo Processed dataset.
  141. */
  142. void AlignWithAkvoDataset( const std::string& file ) {
  143. AlignWithAkvoDataset( YAML::LoadFile(file) );
  144. }
  145. /**
  146. * Assign transmiter coils
  147. */
  148. inline void PushCoil( const std::string& label, std::shared_ptr<PolygonalWireAntenna> ant ) {
  149. TxRx[label] = ant;
  150. }
  151. // ==================== INQUIRY =======================
  152. /**
  153. * @return std::shared_ptr<LayeredEarthEM>
  154. */
  155. inline std::shared_ptr<LayeredEarthEM> GetSigmaModel ( ) {
  156. return SigmaModel;
  157. } // ----- end of method KernelV0::get_SigmaModel -----
  158. /**
  159. * @return the kernel matrix
  160. */
  161. inline MatrixXcr GetKernel ( ) {
  162. return Kern;
  163. }
  164. /**
  165. * @return the integration tolerance
  166. */
  167. inline Real GetTolerance ( ) {
  168. return tol;
  169. }
  170. /**
  171. * @return the layer interfaces
  172. */
  173. inline VectorXr GetInterfaces ( ) {
  174. return Interfaces;
  175. }
  176. /**
  177. * @return the pulse peak current
  178. */
  179. inline VectorXr GetPulseCurrent ( ) {
  180. return PulseI;
  181. }
  182. /**
  183. * @param[in] value the 1D-EM model used for calculations
  184. */
  185. inline void SetLayeredEarthEM ( std::shared_ptr< LayeredEarthEM > value ) {
  186. SigmaModel = value;
  187. return ;
  188. } // ----- end of method KernelV0::set_SigmaModel -----
  189. /**
  190. * @param[in] size the size of the volume to be integrated
  191. */
  192. inline void SetIntegrationSize ( const Vector3r& size ) {
  193. Size = size;
  194. return ;
  195. } // ----- end of method KernelV0::SetIntegrationSize -----
  196. /**
  197. * @param[in] type The type of Hankel transform that will be used.
  198. */
  199. inline void SetHankelTransformType ( const HANKELTRANSFORMTYPE& type ) {
  200. HankelType = type;
  201. return ;
  202. } // ----- end of method KernelV0::SetIntegrationOrigin -----
  203. /**
  204. * @param[in] min is the minimum leaf level, defaults to 0
  205. */
  206. inline void SetMinLevel ( const int& min ) {
  207. minLevel = min;
  208. return ;
  209. } // ----- end of method KernelV0::SetMinLevel -----
  210. /**
  211. * @param[in] max is the maximum leaf level, defaults to 12
  212. */
  213. inline void SetMaxLevel ( const int& max ) {
  214. maxLevel = max;
  215. return ;
  216. } // ----- end of method KernelV0::SetMaxLevel -----
  217. /**
  218. * @param[in] origin The origin location (corner) for the integration volume
  219. */
  220. inline void SetIntegrationOrigin ( const Vector3r& origin ) {
  221. Origin = origin;
  222. return ;
  223. } // ----- end of method KernelV0::SetIntegrationOrigin -----
  224. /**
  225. * @param[in] Amps is the current for each pulse moment
  226. */
  227. inline void SetPulseCurrent ( const VectorXr& Amps ) {
  228. PulseI = Amps;
  229. return ;
  230. } // ----- end of method KernelV0::SetIntegrationOrigin -----
  231. /**
  232. * Sets the temperature, which has implications in calculation of \f$ M_N^{(0)}\f$. Units in
  233. * Kelvin.
  234. */
  235. inline void SetTemperature(const Real& tempK) {
  236. Temperature = tempK;
  237. }
  238. /**
  239. * Sets the tolerance to use for making the adaptive mesh
  240. * @param[in] ttol is the tolerance to use
  241. */
  242. inline void SetTolerance(const Real& ttol) {
  243. tol = ttol;
  244. }
  245. /**
  246. * @param[in] taup sets the pulse duration
  247. */
  248. inline void SetPulseDuration(const Real& taup) {
  249. Taup = taup;
  250. }
  251. inline Real GetPulseDuration( ) {
  252. return Taup;
  253. }
  254. inline void SetDepthLayerInterfaces( const VectorXr& iface ){
  255. Interfaces = iface;
  256. }
  257. /**
  258. * Returns the name of the underlying class, similiar to Python's type
  259. * @return string of class name
  260. */
  261. virtual inline std::string GetName() const {
  262. return CName;
  263. }
  264. protected:
  265. // ==================== LIFECYCLE =======================
  266. /** Copy is disabled */
  267. KernelV0( const KernelV0& ) = delete;
  268. private:
  269. /**
  270. * Returns the kernel value for an input prism
  271. */
  272. VectorXcr f( const Vector3r& r, const Real& volume , const Vector3cr& Ht, const Vector3cr& Hr);
  273. // Complex ComputeV0Cell(const EllipticB& EBT, const EllipticB& EBR,
  274. // const Real& sintheta, const Real& phase, const Real& Mn0Abs,
  275. // const Real& vol);
  276. EllipticB EllipticFieldRep (const Vector3cr& B, const Vector3r& B0hat);
  277. Vector3r ComputeMn0(const Real& Porosity, const Vector3r& B0);
  278. void IntegrateOnOctreeGrid( bool vtkOutput=false );
  279. /**
  280. * Recursive call to integrate a function on an adaptive Octree Grid.
  281. * For efficiency's sake the octree grid is not stored, as only the
  282. * integral (sum) is of interest. The logic for grid refinement is based
  283. * on an Octree representation of the domain. If an Octree representation
  284. * of the kernel is desired, call alternative version @see EvaluateKids2
  285. * @param[in] size gives the domain size, in metres
  286. * @param[in] level gives the current level of the octree grid, call with 0 initially
  287. * @param[in] cpos is the centre position of the parent cuboid
  288. */
  289. void EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
  290. const VectorXcr& parentVal );
  291. #ifdef LEMMAUSEVTK
  292. /**
  293. * Same functionality as @see EvaluateKids, but includes generation of a VTK
  294. * HyperOctree, which is useful for visualization.
  295. */
  296. void EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
  297. const VectorXcr& parentVal, vtkHyperTreeGrid* octree, vtkHyperTreeCursor* curse );
  298. void GetPosition( vtkHyperTreeCursor* Cursor, Real* p );
  299. #endif
  300. // ==================== DATA MEMBERS =========================
  301. int ilay;
  302. int nleaves;
  303. int minLevel=0;
  304. int maxLevel=12;
  305. Real VOLSUM;
  306. Real tol=1e-11;
  307. Real Temperature=283.;
  308. Real Taup = .020; // Sec
  309. Real Larmor;
  310. Vector3r Size;
  311. Vector3r Origin;
  312. VectorXr PulseI;
  313. VectorXr Interfaces;
  314. MatrixXcr Kern;
  315. HANKELTRANSFORMTYPE HankelType=ANDERSON801;
  316. std::shared_ptr< LayeredEarthEM > SigmaModel = nullptr;
  317. std::shared_ptr< FieldPoints > cpoints = nullptr;
  318. std::map< std::string , std::shared_ptr< PolygonalWireAntenna > > TxRx;
  319. std::map< std::string , std::shared_ptr< EMEarth1D > > EMEarths;
  320. #ifdef LEMMAUSEVTK
  321. std::map< int, VectorXcr > LeafDict; // kernel sum for each q
  322. std::map< int, VectorXcr > LeafHt; // Transmitter field
  323. std::map< int, VectorXcr > LeafHr; // Receiver field
  324. std::map< int, int > LeafDictIdx; // index
  325. std::map< int, Real > LeafDictErr; // error value
  326. #endif
  327. #ifdef HAVE_BOOST_PROGRESS
  328. boost::progress_display* disp;
  329. int percent_done;
  330. #endif
  331. // Physical constants and conversion factors
  332. static constexpr Real GAMMA = 2.67518e8; // MKS units
  333. static constexpr Real INVSQRT2 = 0.70710678118654746; // 1/sqrt(2)
  334. static constexpr Real HBAR = 1.05457148e-34; // m2 kg / s
  335. static constexpr Real NH2O = 6.692e28; // [m^3]
  336. static constexpr Real KB = 1.3805e-23; // m^2 kg s-2 K-1
  337. static constexpr Real CHI_N = 3.29e-3; // MKS units
  338. /** ASCII string representation of the class name */
  339. static constexpr auto CName = "KernelV0";
  340. }; // ----- end of class KernelV0 -----
  341. } // ----- end of namespace Lemma ----
  342. /* vim: set tabstop=4 expandtab */
  343. /* vim: set filetype=cpp */
  344. #endif // ----- #ifndef KERNELV0_INC -----