Surface NMR forward modelling
Ви не можете вибрати більше 25 тем Теми мають розпочинатися з літери або цифри, можуть містити дефіси (-) і не повинні перевищувати 35 символів.

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