Surface NMR forward modelling
您最多选择25个主题 主题必须以字母或数字开头,可以包含连字符 (-),并且长度不得超过35个字符

KernelV0.h 15KB

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