c113545f7e7bed67ac69171bc3c2eea2c0a26086
[Merlin.git] / src / KernelV0.cpp
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:25 PM
13 * @version $Id$
14 * @author Trevor Irons (ti)
15 * @email tirons@egi.utah.edu
16 * @copyright Copyright (c) 2016, University of Utah
17 * @copyright Copyright (c) 2016, Lemma Software, LLC
18 */
19
20
21 #include "KernelV0.h"
22 #include "FieldPoints.h"
23
24 namespace Lemma {
25
26 // ==================== FRIEND METHODS =====================
27
28 std::ostream &operator << (std::ostream &stream, const KernelV0 &ob) {
29 stream << ob.Serialize() << "\n---\n"; // End of doc ---
30 return stream;
31 }
32
33 // ==================== LIFECYCLE =======================
34
35 //--------------------------------------------------------------------------------------
36 // Class: KernelV0
37 // Method: KernelV0
38 // Description: constructor (locked)
39 //--------------------------------------------------------------------------------------
40 KernelV0::KernelV0 (const ctor_key&) : LemmaObject( ) {
41
42 } // ----- end of method KernelV0::KernelV0 (constructor) -----
43
44 //--------------------------------------------------------------------------------------
45 // Class: KernelV0
46 // Method: KernelV0
47 // Description: DeSerializing constructor (locked)
48 //--------------------------------------------------------------------------------------
49 KernelV0::KernelV0 (const YAML::Node& node, const ctor_key&) : LemmaObject(node) {
50
51 } // ----- end of method KernelV0::KernelV0 (constructor) -----
52
53 //--------------------------------------------------------------------------------------
54 // Class: KernelV0
55 // Method: NewSP()
56 // Description: public constructor returing a shared_ptr
57 //--------------------------------------------------------------------------------------
58 std::shared_ptr< KernelV0 > KernelV0::NewSP() {
59 return std::make_shared< KernelV0 >( ctor_key() );
60 }
61
62 //--------------------------------------------------------------------------------------
63 // Class: KernelV0
64 // Method: ~KernelV0
65 // Description: destructor (protected)
66 //--------------------------------------------------------------------------------------
67 KernelV0::~KernelV0 () {
68
69 } // ----- end of method KernelV0::~KernelV0 (destructor) -----
70
71 //--------------------------------------------------------------------------------------
72 // Class: KernelV0
73 // Method: Serialize
74 //--------------------------------------------------------------------------------------
75 YAML::Node KernelV0::Serialize ( ) const {
76 YAML::Node node = LemmaObject::Serialize();
77 node.SetTag( GetName() );
78
79 // Coils Transmitters & Receivers
80 for ( auto txm : TxRx) {
81 node[txm.first] = txm.second->Serialize();
82 }
83
84 // LayeredEarthEM
85 node["SigmaModel"] = SigmaModel->Serialize();
86
87 return node;
88 } // ----- end of method KernelV0::Serialize -----
89
90 //--------------------------------------------------------------------------------------
91 // Class: KernelV0
92 // Method: DeSerialize
93 //--------------------------------------------------------------------------------------
94 std::shared_ptr<KernelV0> KernelV0::DeSerialize ( const YAML::Node& node ) {
95 if (node.Tag() != "KernelV0" ) {
96 throw DeSerializeTypeMismatch( "KernelV0", node.Tag());
97 }
98 return std::make_shared< KernelV0 > ( node, ctor_key() );
99 } // ----- end of method KernelV0::DeSerialize -----
100
101 //--------------------------------------------------------------------------------------
102 // Class: KernelV0
103 // Method: DeSerialize
104 //--------------------------------------------------------------------------------------
105 void KernelV0::CalculateK0 (const std::vector< std::string>& Tx, const std::vector<std::string >& Rx,
106 bool vtkOutput ) {
107
108 // All EM calculations will share same field points
109 cpoints = FieldPoints::NewSP();
110 cpoints->SetNumberOfPoints(8);
111 for (auto tx : Tx) {
112 // Set up EMEarth
113 EMEarths[tx] = EMEarth1D::NewSP();
114 EMEarths[tx]->AttachWireAntenna(TxRx[tx]);
115 EMEarths[tx]->AttachLayeredEarthEM(SigmaModel);
116 EMEarths[tx]->AttachFieldPoints( cpoints );
117 EMEarths[tx]->SetFieldsToCalculate(H);
118 // TODO query for method, altough with flat antennae, this is fastest
119 EMEarths[tx]->SetHankelTransformMethod(ANDERSON801);
120 EMEarths[tx]->SetTxRxMode(TX);
121 }
122 for (auto rx : Rx) {
123 if (EMEarths.count(rx)) {
124 EMEarths[rx]->SetTxRxMode(TXRX);
125 } else {
126 EMEarths[rx] = EMEarth1D::NewSP();
127 EMEarths[rx]->AttachWireAntenna(TxRx[rx]);
128 EMEarths[rx]->AttachLayeredEarthEM(SigmaModel);
129 EMEarths[rx]->AttachFieldPoints( cpoints );
130 EMEarths[rx]->SetFieldsToCalculate(H);
131 // TODO query for method, altough with flat antennae, this is fastest
132 EMEarths[rx]->SetHankelTransformMethod(ANDERSON801);
133 EMEarths[rx]->SetTxRxMode(RX);
134 }
135 }
136 IntegrateOnOctreeGrid( 1e-13, vtkOutput );
137 }
138
139 //--------------------------------------------------------------------------------------
140 // Class: KernelV0
141 // Method: IntegrateOnOctreeGrid
142 //--------------------------------------------------------------------------------------
143 void KernelV0::IntegrateOnOctreeGrid( const Real& tolerance, bool vtkOutput) {
144
145 this->tol = tolerance;
146 //Vector3r Size;
147 Size << 200,200,2;
148 //Vector3r Origin;
149 Origin << 0,0,5.0;
150 Vector3r cpos; // centre position
151 //cpos << 100,100,50;
152 cpos = (Size-Origin).array() / 2.;
153 int maxlevel;
154
155 SUM = 0;
156 VOLSUM = 0;
157 nleaves = 0;
158 if (!vtkOutput) {
159 EvaluateKids( Size, 0, cpos, 1e6 );
160 } else {
161 #ifdef LEMMAUSEVTK
162 vtkHyperOctree* oct = vtkHyperOctree::New();
163 oct->SetDimension(3);
164 oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
165 oct->SetSize( Size(0), Size(1), Size(2) );
166 vtkHyperOctreeCursor* curse = oct->NewCellCursor();
167 curse->ToRoot();
168 EvaluateKids2( Size, 0, cpos, 1e6, oct, curse );
169
170 // Fill in leaf data
171 vtkDoubleArray* kr = vtkDoubleArray::New();
172 kr->SetNumberOfComponents(1);
173 kr->SetName("Re($K_0$)");
174 kr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
175 vtkDoubleArray* ki = vtkDoubleArray::New();
176 ki->SetNumberOfComponents(1);
177 ki->SetName("Im($K_0$)");
178 ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
179 vtkDoubleArray* km = vtkDoubleArray::New();
180 km->SetNumberOfComponents(1);
181 km->SetName("mod($K_0$)");
182 km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
183 vtkIntArray* kid = vtkIntArray::New();
184 kid->SetNumberOfComponents(1);
185 kid->SetName("ID");
186 kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
187 //vtkDoubleArray* kerr = vtkDoubleArray::New();
188 // kerr->SetNumberOfComponents(1);
189 // kerr->SetName("Error");
190
191 for (auto leaf : LeafDict) {
192 kr->InsertTuple1( leaf.first, std::real(leaf.second) );
193 ki->InsertTuple1( leaf.first, std::imag(leaf.second) );
194 kid->InsertTuple1( leaf.first, leaf.first );
195 }
196
197 //for (auto leaf : LeafDictErr) {
198 // kerr->InsertTuple1( leaf.first, std::imag(leaf.second) );
199 //}
200
201 oct->GetLeafData()->AddArray(kr);
202 oct->GetLeafData()->AddArray(ki);
203 oct->GetLeafData()->AddArray(km);
204 oct->GetLeafData()->AddArray(kid);
205 //oct->GetLeafData()->AddArray(kerr);
206
207 auto write = vtkXMLHyperOctreeWriter::New();
208 //write.SetDataModeToAscii()
209 write->SetInputData(oct);
210 write->SetFileName("octree.vto");
211 write->Write();
212 write->Delete();
213
214 //kerr->Delete();
215 kid->Delete();
216 kr->Delete();
217 ki->Delete();
218 km->Delete();
219 curse->Delete();
220 oct->Delete();
221 #else
222 throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
223 #endif
224
225 }
226 std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" << Size(0)*Size(1)*Size(2)
227 << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) << std::endl;
228 std::cout << "nleaves\t" << nleaves << std::endl;
229 std::cout << "KSUM\t" << SUM << std::endl;
230
231 }
232
233 //--------------------------------------------------------------------------------------
234 // Class: KernelV0
235 // Method: f
236 //--------------------------------------------------------------------------------------
237 Complex KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
238 //return Complex(volume*Ht.dot(Hr));
239 return ComputeV0Cell(MU0*Ht, MU0*Hr, volume, 1.0);
240 }
241
242 //--------------------------------------------------------------------------------------
243 // Class: KernelV0
244 // Method: ComputeV0Cell
245 //--------------------------------------------------------------------------------------
246 Complex KernelV0::ComputeV0Cell(const Vector3cr& Bt,
247 const Vector3cr& Br, const Real& vol, const Real& phi) {
248
249 // Compute the elliptic fields
250 Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
251 Vector3r B0 = SigmaModel->GetMagneticField();
252
253 // Elliptic representation
254 EllipticB EBT = EllipticFieldRep(Bt, B0hat);
255 EllipticB EBR = EllipticFieldRep(Br, B0hat);
256
257 // Compute Mn0
258 Vector3r Mn0 = ComputeMn0(phi, B0);
259 Real Mn0Abs = Mn0.norm();
260
261 Real Taup = 0.020; // s
262 Real Ip = 10; // A
263
264 // Compute the tipping angle
265 Real sintheta = std::sin(0.5*GAMMA*Ip*Taup*std::abs(EBT.alpha-EBT.beta));
266
267 // Compute phase delay, TODO add transmiiter phase and delay time phase!
268 Real phase = EBR.zeta+EBT.zeta;
269
270 return ComputeV0Cell(EBT, EBR, sintheta, phase, Mn0Abs, vol);
271 }
272
273 //--------------------------------------------------------------------------------------
274 // Class: KernelV0
275 // Method: ComputeV0Cell
276 //--------------------------------------------------------------------------------------
277 Complex KernelV0::ComputeV0Cell(const EllipticB& EBT, const EllipticB& EBR,
278 const Real& sintheta, const Real& phase, const Real& Mn0Abs,
279 const Real& vol) {
280
281 Vector3r B0hat = {1,0,0};
282 Real Larmor = 2500.*2.*PI;
283
284 // earth response of receiver adjoint field
285 Complex ejztr = std::exp(Complex(0, EBR.zeta + EBT.zeta));
286
287 Complex PhaseTerm = EBR.bhat.dot(EBT.bhat) +
288 (B0hat.dot(EBR.bhat.cross(EBT.bhat) ));
289 return -vol*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
290 }
291
292 //--------------------------------------------------------------------------------------
293 // Class: KernelV0
294 // Method: ComputeV0Cell
295 //--------------------------------------------------------------------------------------
296 Vector3r KernelV0::ComputeMn0(const Real& Porosity, const Vector3r& B0) {
297 Real chi_n = NH2O*((GAMMA*GAMMA*HBAR*HBAR)/(4.*KB*Temperature));
298 return chi_n*Porosity*B0;
299 }
300
301 //--------------------------------------------------------------------------------------
302 // Class: KernelV0
303 // Method: ComputeV0Cell
304 //--------------------------------------------------------------------------------------
305 EllipticB KernelV0::EllipticFieldRep (const Vector3cr& B, const Vector3r& B0hat) {
306 EllipticB ElipB = EllipticB();
307 Vector3cr Bperp = B.array() - B0hat.dot(B)*B0hat.array();
308 Real BperpNorm = Bperp.norm();
309 Complex Bp2 = Bperp.transpose() * Bperp;
310 VectorXcr iB0 = Complex(0,1)*B0hat.cast<Complex>().array();
311 ElipB.eizt = std::sqrt(Bp2 / std::abs(Bp2));
312 ElipB.alpha = INVSQRT2*std::sqrt(BperpNorm*BperpNorm + std::abs(Bp2));
313 ElipB.beta = sgn(std::real(iB0.dot(Bperp.cross(Bperp.conjugate())))) *
314 (INVSQRT2)*std::sqrt(std::abs(BperpNorm*BperpNorm-std::abs(Bp2)));
315 ElipB.bhat = ((Real)1./ElipB.alpha)*(((Real)1./ElipB.eizt)*Bperp.array()).real().array();
316 ElipB.bhatp = B0hat.cross(ElipB.bhat);
317 ElipB.zeta = std::real(std::log(ElipB.eizt)/Complex(0,1));
318 return ElipB;
319 }
320
321 //--------------------------------------------------------------------------------------
322 // Class: KernelV0
323 // Method: EvaluateKids
324 //--------------------------------------------------------------------------------------
325 bool KernelV0::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
326 const Complex& parentVal ) {
327
328 std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
329 std::cout.flush();
330
331 // Next level step, interested in one level below
332 // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
333 Vector3r step = size.array() / (Real)(1 << (level+1) );
334 Vector3r step2 = size.array() / (Real)(1 << (level+2) );
335
336 Real vol = (step2(0)*step2(1)*step2(2)); // volume of each child
337
338 Vector3r pos = cpos - step/2.;
339 Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
340 0, 0, 0,
341 step[0], 0, 0,
342 0, step[1], 0,
343 step[0], step[1], 0,
344 0, 0, step[2],
345 step[0], 0, step[2],
346 0, step[1], step[2],
347 step[0], step[1], step[2] ).finished();
348
349 VectorXcr kvals(8); // individual kernel vals
350 cpoints->ClearFields();
351 for (int ichild=0; ichild<8; ++ichild) {
352 Vector3r cp = pos; // Eigen complains about combining these
353 cp += posadd.row(ichild);
354 cpoints->SetLocation( ichild, cp );
355 }
356
357 Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
358 Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
359 //Eigen::Matrix< Complex, 8, 3 > Bt;
360 for ( auto EMCalc : EMEarths ) {
361 //EMCalc->GetFieldPoints()->ClearFields();
362 EMCalc.second->CalculateWireAntennaFields();
363 switch (EMCalc.second->GetTxRxMode()) {
364 case TX:
365 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
366 break;
367 case RX:
368 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
369 break;
370 case TXRX:
371 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
372 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
373 break;
374 default:
375 break;
376 }
377 }
378
379 for (int ichild=0; ichild<8; ++ichild) {
380 Vector3r cp = pos; // Eigen complains about combining these
381 cp += posadd.row(ichild);
382 kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
383 }
384
385 Complex ksum = kvals.sum(); // Kernel sum
386 // Evaluate whether or not furthur splitting is needed
387 if ( std::abs(ksum - parentVal) > tol || level < 2 ) {
388 for (int ichild=0; ichild<8; ++ichild) {
389 Vector3r cp = pos; // Eigen complains about combining these
390 cp += posadd.row(ichild);
391 bool isleaf = EvaluateKids( size, level+1, cp, kvals(ichild) );
392 if (isleaf) { // Include result in final integral
393 SUM += ksum;
394 VOLSUM += 8.*vol;
395 nleaves += 1;
396 }
397 }
398 return false; // not leaf
399 }
400 // Save here instead?
401 return true; // leaf
402 }
403
404 #ifdef LEMMAUSEVTK
405 //--------------------------------------------------------------------------------------
406 // Class: KernelV0
407 // Method: EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
408 //--------------------------------------------------------------------------------------
409 bool KernelV0::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
410 const Complex& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
411
412 std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
413 std::cout.flush();
414
415 // Next level step, interested in one level below
416 // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
417 Vector3r step = size.array() / (Real)(1 << (level+1) );
418 Vector3r step2 = size.array() / (Real)(1 << (level+2) );
419
420 Real vol = (step2(0)*step2(1)*step2(2)); // volume of each child
421
422 Vector3r pos = cpos - step/2.;
423 Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
424 0, 0, 0,
425 step[0], 0, 0,
426 0, step[1], 0,
427 step[0], step[1], 0,
428 0, 0, step[2],
429 step[0], 0, step[2],
430 0, step[1], step[2],
431 step[0], step[1], step[2] ).finished();
432
433 VectorXcr kvals(8); // individual kernel vals
434 cpoints->ClearFields();
435 for (int ichild=0; ichild<8; ++ichild) {
436 Vector3r cp = pos; // Eigen complains about combining these
437 cp += posadd.row(ichild);
438 cpoints->SetLocation( ichild, cp );
439 }
440
441 Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
442 Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
443 for ( auto EMCalc : EMEarths ) {
444 //EMCalc->GetFieldPoints()->ClearFields();
445 EMCalc.second->CalculateWireAntennaFields();
446 switch (EMCalc.second->GetTxRxMode()) {
447 case TX:
448 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
449 break;
450 case RX:
451 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
452 break;
453 case TXRX:
454 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
455 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
456 break;
457 default:
458 break;
459 }
460 }
461
462 for (int ichild=0; ichild<8; ++ichild) {
463 Vector3r cp = pos; // Eigen complains about combining these
464 cp += posadd.row(ichild);
465 kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
466 }
467
468 Complex ksum = kvals.sum(); // Kernel sum
469 // Evaluate whether or not furthur splitting is needed
470 if ( std::abs(ksum - parentVal) > tol || level < 2 ) {
471 oct->SubdivideLeaf(curse);
472 for (int ichild=0; ichild<8; ++ichild) {
473 curse->ToChild(ichild);
474 Vector3r cp = pos; // Eigen complains about combining these
475 cp += posadd.row(ichild);
476 /* Test for position via alternative means
477 Real p[3];
478 GetPosition(curse, p);
479 std::cout << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
480 << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
481 */
482 bool isleaf = EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
483 if (isleaf) { // Include result in final integral
484 LeafDict[curse->GetLeafId()] = kvals(ichild); // VTK
485 SUM += ksum;
486 VOLSUM += 8*vol;
487 nleaves += 1;
488 }
489 curse->ToParent();
490 }
491 return false; // not leaf
492 }
493 return true; // leaf
494 }
495
496 //--------------------------------------------------------------------------------------
497 // Class: KernelV0
498 // Method: GetPosition
499 //--------------------------------------------------------------------------------------
500 void KernelV0::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
501 Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
502 //step = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
503 p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
504 p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
505 p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
506 }
507
508 #endif
509
510 } // ---- end of namespace Lemma ----
511
512 /* vim: set tabstop=4 expandtab */
513 /* vim: set filetype=cpp */
514