Kernel calculation seems to be working.
[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 // Set up
109 Larmor = SigmaModel->GetMagneticFieldMagnitude()*GAMMA; // in rad 2246.*2.*PI;
110
111 // All EM calculations will share same field points
112 cpoints = FieldPoints::NewSP();
113 cpoints->SetNumberOfPoints(8);
114 for (auto tx : Tx) {
115 // Set up EMEarth
116 EMEarths[tx] = EMEarth1D::NewSP();
117 EMEarths[tx]->AttachWireAntenna(TxRx[tx]);
118 EMEarths[tx]->AttachLayeredEarthEM(SigmaModel);
119 EMEarths[tx]->AttachFieldPoints( cpoints );
120 EMEarths[tx]->SetFieldsToCalculate(H);
121 // TODO query for method, altough with flat antennae, this is fastest
122 EMEarths[tx]->SetHankelTransformMethod(ANDERSON801);
123 EMEarths[tx]->SetTxRxMode(TX);
124 }
125 for (auto rx : Rx) {
126 if (EMEarths.count(rx)) {
127 EMEarths[rx]->SetTxRxMode(TXRX);
128 } else {
129 EMEarths[rx] = EMEarth1D::NewSP();
130 EMEarths[rx]->AttachWireAntenna(TxRx[rx]);
131 EMEarths[rx]->AttachLayeredEarthEM(SigmaModel);
132 EMEarths[rx]->AttachFieldPoints( cpoints );
133 EMEarths[rx]->SetFieldsToCalculate(H);
134 // TODO query for method, altough with flat antennae, this is fastest
135 EMEarths[rx]->SetHankelTransformMethod(ANDERSON801);
136 EMEarths[rx]->SetTxRxMode(RX);
137 }
138 }
139
140 std::cout << "Calculating K0 kernel\n";
141 MatrixXcr Kern = MatrixXcr::Zero( Interfaces.size() - 1, PulseI.size() );
142 for (int ilay=0; ilay< Interfaces.size()-1; ++ilay) {
143 for (int iq=0; iq< PulseI.size()-1; ++iq) {
144 std::cout << "Layer " << ilay << " q " << iq << std::endl;
145 Size(2) = Interfaces(ilay+1) - Interfaces(ilay);
146 Origin(2) = Interfaces(ilay);
147 Ip = PulseI(iq);
148 Kern(ilay, iq) = IntegrateOnOctreeGrid( ilay, iq, vtkOutput );
149 }
150 }
151 std::cout << "\rFinished KERNEL\n";
152 std::cout << Kern << std::endl;
153 //IntegrateOnOctreeGrid( vtkOutput );
154 }
155
156 //--------------------------------------------------------------------------------------
157 // Class: KernelV0
158 // Method: IntegrateOnOctreeGrid
159 //--------------------------------------------------------------------------------------
160 Complex KernelV0::IntegrateOnOctreeGrid( const int& ilay, const int& iq, bool vtkOutput) {
161
162 Vector3r cpos = (Size-Origin).array() / 2.;
163
164 SUM = 0;
165 VOLSUM = 0;
166 nleaves = 0;
167 if (!vtkOutput) {
168 EvaluateKids( Size, 0, cpos, 1e6 );
169 } else {
170 #ifdef LEMMAUSEVTK
171 vtkHyperOctree* oct = vtkHyperOctree::New();
172 oct->SetDimension(3);
173 oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
174 oct->SetSize( Size(0), Size(1), Size(2) );
175 vtkHyperOctreeCursor* curse = oct->NewCellCursor();
176 curse->ToRoot();
177 EvaluateKids2( Size, 0, cpos, 1e6, oct, curse );
178
179 // Fill in leaf data
180 vtkDoubleArray* kr = vtkDoubleArray::New();
181 kr->SetNumberOfComponents(1);
182 kr->SetName("Re($K_0$)");
183 kr->SetNumberOfTuples( oct->GetNumberOfLeaves() );
184 vtkDoubleArray* ki = vtkDoubleArray::New();
185 ki->SetNumberOfComponents(1);
186 ki->SetName("Im($K_0$)");
187 ki->SetNumberOfTuples( oct->GetNumberOfLeaves() );
188 vtkDoubleArray* km = vtkDoubleArray::New();
189 km->SetNumberOfComponents(1);
190 km->SetName("mod($K_0$)");
191 km->SetNumberOfTuples( oct->GetNumberOfLeaves() );
192 vtkIntArray* kid = vtkIntArray::New();
193 kid->SetNumberOfComponents(1);
194 kid->SetName("ID");
195 kid->SetNumberOfTuples( oct->GetNumberOfLeaves() );
196 vtkIntArray* kerr = vtkIntArray::New();
197 kerr->SetNumberOfComponents(1);
198 kerr->SetName("nleaf");
199
200 for (auto leaf : LeafDict) {
201 kr->InsertTuple1( leaf.first, std::real(leaf.second) );
202 ki->InsertTuple1( leaf.first, std::imag(leaf.second) );
203 km->InsertTuple1( leaf.first, std::abs(leaf.second) );
204 kid->InsertTuple1( leaf.first, leaf.first );
205 }
206
207 for (auto leaf : LeafDictIdx) {
208 kerr->InsertTuple1( leaf.first, leaf.second );
209 }
210
211 oct->GetLeafData()->AddArray(kr);
212 oct->GetLeafData()->AddArray(ki);
213 oct->GetLeafData()->AddArray(km);
214 oct->GetLeafData()->AddArray(kid);
215 oct->GetLeafData()->AddArray(kerr);
216
217 auto write = vtkXMLHyperOctreeWriter::New();
218 //write.SetDataModeToAscii()
219 write->SetInputData(oct);
220 std::string fname = std::string("octree-") + to_string(ilay)
221 + std::string("-") + to_string(iq) + std::string(".vto");
222 write->SetFileName(fname.c_str());
223 write->Write();
224 write->Delete();
225
226 //kerr->Delete();
227 kid->Delete();
228 kr->Delete();
229 ki->Delete();
230 km->Delete();
231 curse->Delete();
232 oct->Delete();
233 #else
234 throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
235 #endif
236
237 }
238 std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" << Size(0)*Size(1)*Size(2)
239 << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) << std::endl;
240 std::cout << "nleaves\t" << nleaves << std::endl;
241 std::cout << "KSUM\t" << SUM << std::endl;
242
243 return SUM;
244 }
245
246 //--------------------------------------------------------------------------------------
247 // Class: KernelV0
248 // Method: f
249 //--------------------------------------------------------------------------------------
250 Complex KernelV0::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
251 //return Complex(volume*Ht.dot(Hr));
252 return ComputeV0Cell(MU0*Ht, MU0*Hr, volume, 1.0);
253 }
254
255 //--------------------------------------------------------------------------------------
256 // Class: KernelV0
257 // Method: ComputeV0Cell
258 //--------------------------------------------------------------------------------------
259 Complex KernelV0::ComputeV0Cell(const Vector3cr& Bt,
260 const Vector3cr& Br, const Real& vol, const Real& phi) {
261
262 // Compute the elliptic fields
263 Vector3r B0hat = SigmaModel->GetMagneticFieldUnitVector();
264 Vector3r B0 = SigmaModel->GetMagneticField();
265
266 // Elliptic representation
267 EllipticB EBT = EllipticFieldRep(Bt, B0hat);
268 EllipticB EBR = EllipticFieldRep(Br, B0hat);
269
270 // Compute Mn0
271 Vector3r Mn0 = ComputeMn0(phi, B0);
272 Real Mn0Abs = Mn0.norm();
273
274 // Compute the tipping angle
275 Real sintheta = std::sin(0.5*GAMMA*Ip*Taup*std::abs(EBT.alpha-EBT.beta));
276
277 // Compute phase delay, TODO add transmiiter phase and delay time phase!
278 Real phase = EBR.zeta+EBT.zeta;
279
280 return ComputeV0Cell(EBT, EBR, sintheta, phase, Mn0Abs, vol);
281 }
282
283 //--------------------------------------------------------------------------------------
284 // Class: KernelV0
285 // Method: ComputeV0Cell
286 //--------------------------------------------------------------------------------------
287 Complex KernelV0::ComputeV0Cell(const EllipticB& EBT, const EllipticB& EBR,
288 const Real& sintheta, const Real& phase, const Real& Mn0Abs,
289 const Real& vol) {
290
291 Vector3r B0hat = {1,0,0};
292
293 // earth response of receiver adjoint field
294 Complex ejztr = std::exp(Complex(0, EBR.zeta + EBT.zeta));
295
296 Complex PhaseTerm = EBR.bhat.dot(EBT.bhat) +
297 (B0hat.dot(EBR.bhat.cross(EBT.bhat) ));
298 return -vol*Complex(0,Larmor)*Mn0Abs*(EBR.alpha+EBR.beta)*ejztr*sintheta*PhaseTerm;
299 }
300
301 //--------------------------------------------------------------------------------------
302 // Class: KernelV0
303 // Method: ComputeV0Cell
304 //--------------------------------------------------------------------------------------
305 Vector3r KernelV0::ComputeMn0(const Real& Porosity, const Vector3r& B0) {
306 Real chi_n = NH2O*((GAMMA*GAMMA*HBAR*HBAR)/(4.*KB*Temperature));
307 return chi_n*Porosity*B0;
308 }
309
310 //--------------------------------------------------------------------------------------
311 // Class: KernelV0
312 // Method: ComputeV0Cell
313 //--------------------------------------------------------------------------------------
314 EllipticB KernelV0::EllipticFieldRep (const Vector3cr& B, const Vector3r& B0hat) {
315 EllipticB ElipB = EllipticB();
316 Vector3cr Bperp = B.array() - B0hat.dot(B)*B0hat.array();
317 Real BperpNorm = Bperp.norm();
318 Complex Bp2 = Bperp.transpose() * Bperp;
319 VectorXcr iB0 = Complex(0,1)*B0hat.cast<Complex>().array();
320 ElipB.eizt = std::sqrt(Bp2 / std::abs(Bp2));
321 ElipB.alpha = INVSQRT2*std::sqrt(BperpNorm*BperpNorm + std::abs(Bp2));
322 ElipB.beta = sgn(std::real(iB0.dot(Bperp.cross(Bperp.conjugate())))) *
323 (INVSQRT2)*std::sqrt(std::abs(BperpNorm*BperpNorm-std::abs(Bp2)));
324 ElipB.bhat = ((Real)1./ElipB.alpha)*(((Real)1./ElipB.eizt)*Bperp.array()).real().array();
325 ElipB.bhatp = B0hat.cross(ElipB.bhat);
326 ElipB.zeta = std::real(std::log(ElipB.eizt)/Complex(0,1));
327 return ElipB;
328 }
329
330 //--------------------------------------------------------------------------------------
331 // Class: KernelV0
332 // Method: EvaluateKids
333 //--------------------------------------------------------------------------------------
334 bool KernelV0::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
335 const Complex& parentVal ) {
336
337 std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
338 std::cout.flush();
339
340 // Next level step, interested in one level below
341 // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
342 Vector3r step = size.array() / (Real)(1 << (level+1) );
343 Vector3r step2 = size.array() / (Real)(1 << (level+2) );
344
345 Real vol = (step2(0)*step2(1)*step2(2)); // volume of each child
346
347 Vector3r pos = cpos - step/2.;
348 Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
349 0, 0, 0,
350 step[0], 0, 0,
351 0, step[1], 0,
352 step[0], step[1], 0,
353 0, 0, step[2],
354 step[0], 0, step[2],
355 0, step[1], step[2],
356 step[0], step[1], step[2] ).finished();
357
358 VectorXcr kvals(8); // individual kernel vals
359 cpoints->ClearFields();
360 for (int ichild=0; ichild<8; ++ichild) {
361 Vector3r cp = pos; // Eigen complains about combining these
362 cp += posadd.row(ichild);
363 cpoints->SetLocation( ichild, cp );
364 }
365
366 Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
367 Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
368 //Eigen::Matrix< Complex, 8, 3 > Bt;
369 for ( auto EMCalc : EMEarths ) {
370 //EMCalc->GetFieldPoints()->ClearFields();
371 EMCalc.second->CalculateWireAntennaFields();
372 switch (EMCalc.second->GetTxRxMode()) {
373 case TX:
374 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
375 break;
376 case RX:
377 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
378 break;
379 case TXRX:
380 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
381 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
382 break;
383 default:
384 break;
385 }
386 }
387
388 for (int ichild=0; ichild<8; ++ichild) {
389 Vector3r cp = pos; // Eigen complains about combining these
390 cp += posadd.row(ichild);
391 kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
392 }
393
394 Complex ksum = kvals.sum(); // Kernel sum
395 // Evaluate whether or not furthur splitting is needed
396 if ( std::abs(ksum - parentVal) > tol || level < minLevel && level < maxLevel ) {
397 for (int ichild=0; ichild<8; ++ichild) {
398 Vector3r cp = pos; // Eigen complains about combining these
399 cp += posadd.row(ichild);
400 bool isleaf = EvaluateKids( size, level+1, cp, kvals(ichild) );
401 if (isleaf) { // Include result in final integral
402 SUM += ksum;
403 VOLSUM += 8.*vol;
404 nleaves += 1;
405 }
406 }
407 return false; // not leaf
408 }
409 // Save here instead?
410 return true; // leaf
411 }
412
413 #ifdef LEMMAUSEVTK
414 //--------------------------------------------------------------------------------------
415 // Class: KernelV0
416 // Method: EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
417 //--------------------------------------------------------------------------------------
418 bool KernelV0::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
419 const Complex& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
420
421 std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
422 std::cout.flush();
423
424 // Next level step, interested in one level below
425 // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
426 Vector3r step = size.array() / (Real)(1 << (level+1) );
427 Vector3r step2 = size.array() / (Real)(1 << (level+2) );
428
429 Real vol = (step2(0)*step2(1)*step2(2)); // volume of each child
430
431 Vector3r pos = cpos - step/2.;
432 Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
433 0, 0, 0,
434 step[0], 0, 0,
435 0, step[1], 0,
436 step[0], step[1], 0,
437 0, 0, step[2],
438 step[0], 0, step[2],
439 0, step[1], step[2],
440 step[0], step[1], step[2] ).finished();
441
442 VectorXcr kvals(8); // individual kernel vals
443 cpoints->ClearFields();
444 for (int ichild=0; ichild<8; ++ichild) {
445 Vector3r cp = pos; // Eigen complains about combining these
446 cp += posadd.row(ichild);
447 cpoints->SetLocation( ichild, cp );
448 }
449
450 Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
451 Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
452 for ( auto EMCalc : EMEarths ) {
453 //EMCalc->GetFieldPoints()->ClearFields();
454 EMCalc.second->CalculateWireAntennaFields();
455 switch (EMCalc.second->GetTxRxMode()) {
456 case TX:
457 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
458 break;
459 case RX:
460 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
461 break;
462 case TXRX:
463 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
464 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
465 break;
466 default:
467 break;
468 }
469 }
470
471 for (int ichild=0; ichild<8; ++ichild) {
472 Vector3r cp = pos; // Eigen complains about combining these
473 cp += posadd.row(ichild);
474 kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
475 }
476
477 Complex ksum = kvals.sum(); // Kernel sum
478 // Evaluate whether or not furthur splitting is needed
479 Real err = std::abs(ksum - parentVal);
480 if ( std::abs(ksum - parentVal) > tol || level < minLevel && level < maxLevel ) {
481 oct->SubdivideLeaf(curse);
482 for (int ichild=0; ichild<8; ++ichild) {
483 curse->ToChild(ichild);
484 Vector3r cp = pos; // Eigen complains about combining these
485 cp += posadd.row(ichild);
486 /* Test for position via alternative means
487 Real p[3];
488 GetPosition(curse, p);
489 std::cout << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
490 << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
491 */
492 bool isleaf = EvaluateKids2( size, level+1, cp, kvals(ichild), oct, curse );
493 if (isleaf) { // Include result in final integral
494 LeafDict[curse->GetLeafId()] = kvals(ichild); // VTK
495 LeafDictIdx[curse->GetLeafId()] = nleaves; // VTK
496 SUM += ksum;
497 VOLSUM += 8*vol;
498 nleaves += 1;
499 }
500 curse->ToParent();
501 }
502 return false; // not leaf
503 }
504 return true; // leaf
505 }
506
507 //--------------------------------------------------------------------------------------
508 // Class: KernelV0
509 // Method: GetPosition
510 //--------------------------------------------------------------------------------------
511 void KernelV0::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
512 Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
513 //step = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
514 p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
515 p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
516 p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
517 }
518
519 #endif
520
521 } // ---- end of namespace Lemma ----
522
523 /* vim: set tabstop=4 expandtab */
524 /* vim: set filetype=cpp */
525