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