1d84f72e43480340bfc6b911b350f962c351a705
[Merlin.git] / src / Coupling.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 * @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
20
21 #include "Coupling.h"
22 #include "FieldPoints.h"
23
24 namespace Lemma {
25
26 // ==================== FRIEND METHODS =====================
27
28 std::ostream &operator << (std::ostream &stream, const Coupling &ob) {
29 stream << ob.Serialize() << "\n---\n"; // End of doc ---
30 return stream;
31 }
32
33 // ==================== LIFECYCLE =======================
34
35 //--------------------------------------------------------------------------------------
36 // Class: Coupling
37 // Method: Coupling
38 // Description: constructor (locked)
39 //--------------------------------------------------------------------------------------
40 Coupling::Coupling (const ctor_key&) : LemmaObject( ) {
41
42 } // ----- end of method Coupling::Coupling (constructor) -----
43
44 //--------------------------------------------------------------------------------------
45 // Class: Coupling
46 // Method: Coupling
47 // Description: DeSerializing constructor (locked)
48 //--------------------------------------------------------------------------------------
49 Coupling::Coupling (const YAML::Node& node, const ctor_key&) : LemmaObject(node) {
50
51 } // ----- end of method Coupling::Coupling (constructor) -----
52
53 //--------------------------------------------------------------------------------------
54 // Class: Coupling
55 // Method: NewSP()
56 // Description: public constructor returing a shared_ptr
57 //--------------------------------------------------------------------------------------
58 std::shared_ptr< Coupling > Coupling::NewSP() {
59 return std::make_shared< Coupling >( ctor_key() );
60 }
61
62 //--------------------------------------------------------------------------------------
63 // Class: Coupling
64 // Method: ~Coupling
65 // Description: destructor (protected)
66 //--------------------------------------------------------------------------------------
67 Coupling::~Coupling () {
68
69 } // ----- end of method Coupling::~Coupling (destructor) -----
70
71 //--------------------------------------------------------------------------------------
72 // Class: Coupling
73 // Method: Serialize
74 //--------------------------------------------------------------------------------------
75 YAML::Node Coupling::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 // LayeredEarthEM
84 node["SigmaModel"] = SigmaModel->Serialize();
85
86 node["tol"] = tol;
87 node["minLevel"] = minLevel;
88 node["maxLevel"] = maxLevel;
89
90 return node;
91 } // ----- end of method Coupling::Serialize -----
92
93 //--------------------------------------------------------------------------------------
94 // Class: Coupling
95 // Method: DeSerialize
96 //--------------------------------------------------------------------------------------
97 std::shared_ptr<Coupling> Coupling::DeSerialize ( const YAML::Node& node ) {
98 if (node.Tag() != "Coupling" ) {
99 throw DeSerializeTypeMismatch( "Coupling", node.Tag());
100 }
101 return std::make_shared< Coupling > ( node, ctor_key() );
102 } // ----- end of method Coupling::DeSerialize -----
103
104 //--------------------------------------------------------------------------------------
105 // Class: Coupling
106 // Method: DeSerialize
107 //--------------------------------------------------------------------------------------
108 Complex Coupling::Calculate (const std::vector< std::string>& Tx, const std::vector<std::string >& Rx,
109 bool vtkOutput ) {
110 // All EM calculations will share same field points
111 cpoints = FieldPoints::NewSP();
112 cpoints->SetNumberOfPoints(8);
113 for (auto tx : Tx) {
114 // Set up EMEarth
115 EMEarths[tx] = EMEarth1D::NewSP();
116 EMEarths[tx]->AttachWireAntenna(TxRx[tx]);
117 EMEarths[tx]->AttachLayeredEarthEM(SigmaModel);
118 EMEarths[tx]->AttachFieldPoints( cpoints );
119 EMEarths[tx]->SetFieldsToCalculate(H);
120 // TODO query for method, altough with flat antennae, this is fastest
121 EMEarths[tx]->SetHankelTransformMethod(ANDERSON801);
122 EMEarths[tx]->SetTxRxMode(TX);
123 TxRx[tx]->SetCurrent(1.);
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 TxRx[rx]->SetCurrent(1.);
138 }
139 }
140 SUM = 0;
141 IntegrateOnOctreeGrid( vtkOutput );
142 std::cout << "\nFinished KERNEL\n";
143 return SUM;
144 }
145
146 //--------------------------------------------------------------------------------------
147 // Class: Coupling
148 // Method: IntegrateOnOctreeGrid
149 //--------------------------------------------------------------------------------------
150 void Coupling::IntegrateOnOctreeGrid( bool vtkOutput) {
151
152 Vector3r cpos = Origin + Size/2.;
153
154 VOLSUM = 0;
155 nleaves = 0;
156 if (!vtkOutput) {
157 EvaluateKids( Size, 0, cpos, Complex(100.));
158 } else {
159 #ifdef LEMMAUSEVTK
160 vtkHyperOctree* oct = vtkHyperOctree::New();
161 oct->SetDimension(3);
162 oct->SetOrigin( Origin(0), Origin(1), Origin(2) );
163 oct->SetSize( Size(0), Size(1), Size(2) );
164 vtkHyperOctreeCursor* curse = oct->NewCellCursor();
165 curse->ToRoot();
166 EvaluateKids2( Size, 0, cpos, VectorXcr::Ones(PulseI.size()), oct, curse );
167
168 for (int iq=0; iq<PulseI.size(); ++iq) {
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 vtkIntArray* kerr = vtkIntArray::New();
188 kerr->SetNumberOfComponents(1);
189 kerr->SetName("nleaf");
190
191 //Real LeafVol(0);
192 for (auto leaf : LeafDict) {
193 kr->InsertTuple1( leaf.first, std::real(leaf.second(iq)) );
194 ki->InsertTuple1( leaf.first, std::imag(leaf.second(iq)) );
195 km->InsertTuple1( leaf.first, std::abs(leaf.second(iq)) );
196 kid->InsertTuple1( leaf.first, leaf.first );
197 //LeafVol += std::real(leaf.second);
198 }
199 //std::cout << "\n\nLeafVol=" << LeafVol << std::endl;
200
201 for (auto leaf : LeafDictIdx) {
202 kerr->InsertTuple1( leaf.first, leaf.second );
203 }
204
205 auto kri = oct->GetLeafData()->AddArray(kr);
206 auto kii = oct->GetLeafData()->AddArray(ki);
207 auto kmi = oct->GetLeafData()->AddArray(km);
208 auto kidi = oct->GetLeafData()->AddArray(kid);
209 auto keri = oct->GetLeafData()->AddArray(kerr);
210
211 auto write = vtkXMLHyperOctreeWriter::New();
212 //write.SetDataModeToAscii()
213 write->SetInputData(oct);
214 std::string fname = std::string("octree-") + to_string(ilay)
215 + std::string("-") + to_string(iq) + std::string(".vto");
216 write->SetFileName(fname.c_str());
217 write->Write();
218 write->Delete();
219
220 oct->GetLeafData()->RemoveArray( kri );
221 oct->GetLeafData()->RemoveArray( kii );
222 oct->GetLeafData()->RemoveArray( kmi );
223 oct->GetLeafData()->RemoveArray( kidi );
224 oct->GetLeafData()->RemoveArray( keri );
225
226 kerr->Delete();
227 kid->Delete();
228 kr->Delete();
229 ki->Delete();
230 km->Delete();
231
232 }
233
234 curse->Delete();
235 oct->Delete();
236 #else
237 throw std::runtime_error("IntegrateOnOctreeGrid with vtkOutput requires Lemma with VTK support");
238 #endif
239
240 }
241 std::cout << "\nVOLSUM=" << VOLSUM << "\tActual=" << Size(0)*Size(1)*Size(2)
242 << "\tDifference=" << VOLSUM - (Size(0)*Size(1)*Size(2)) << std::endl;
243 }
244
245 //--------------------------------------------------------------------------------------
246 // Class: Coupling
247 // Method: f
248 //--------------------------------------------------------------------------------------
249 Complex Coupling::f( const Vector3r& r, const Real& volume, const Vector3cr& Ht, const Vector3cr& Hr ) {
250 return volume*Ht.dot(Hr);
251 }
252
253 //--------------------------------------------------------------------------------------
254 // Class: Coupling
255 // Method: EvaluateKids
256 //--------------------------------------------------------------------------------------
257 void Coupling::EvaluateKids( const Vector3r& size, const int& level, const Vector3r& cpos,
258 const Complex& parentVal ) {
259
260 std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
261 //std::cout.flush();
262
263 // Next level step, interested in one level below
264 // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
265 Vector3r step = size.array() / (Real)(1 << (level+1) );
266 Real vol = (step(0)*step(1)*step(2)); // volume of each child
267
268 Vector3r pos = cpos - step/2.;
269 Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
270 0, 0, 0,
271 step[0], 0, 0,
272 0, step[1], 0,
273 step[0], step[1], 0,
274 0, 0, step[2],
275 step[0], 0, step[2],
276 0, step[1], step[2],
277 step[0], step[1], step[2] ).finished();
278
279 VectorXcr kvals(8); // individual kernel vals
280 cpoints->ClearFields();
281 for (int ichild=0; ichild<8; ++ichild) {
282 Vector3r cp = pos; // Eigen complains about combining these
283 cp += posadd.row(ichild);
284 cpoints->SetLocation( ichild, cp );
285 }
286
287 Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
288 Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
289 for ( auto EMCalc : EMEarths ) {
290 EMCalc.second->GetFieldPoints()->ClearFields();
291 EMCalc.second->CalculateWireAntennaFields();
292 switch (EMCalc.second->GetTxRxMode()) {
293 case TX:
294 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
295 break;
296 case RX:
297 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
298 break;
299 case TXRX:
300 Ht += EMCalc.second->GetFieldPoints()->GetHfield(0);
301 Hr += EMCalc.second->GetFieldPoints()->GetHfield(0);
302 break;
303 default:
304 break;
305 }
306 }
307
308 for (int ichild=0; ichild<8; ++ichild) {
309 Vector3r cp = pos; // Eigen complains about combining these
310 cp += posadd.row(ichild);
311 kvals(ichild) = f(cp, vol, Ht.col(ichild), Hr.col(ichild));
312 }
313
314 Complex ksum = kvals.sum(); // Kernel sum
315 // Evaluate whether or not furthur splitting is needed
316 if ( std::abs(ksum-parentVal) > tol || level < minLevel && level < maxLevel ) {
317 // Not a leaf dive further in
318 for (int ichild=0; ichild<8; ++ichild) {
319 Vector3r cp = pos; // Eigen complains about combining these
320 cp += posadd.row(ichild);
321 EvaluateKids( size, level+1, cp, kvals(ichild) );
322 }
323 return; // not leaf
324 }
325 // implicit else, is a leaf
326 SUM += ksum;
327 VOLSUM += 8.*vol;
328 nleaves += 1; // could say += 8 just as fairly
329 return; // is leaf
330 }
331
332 #ifdef LEMMAUSEVTK
333 //--------------------------------------------------------------------------------------
334 // Class: Coupling
335 // Method: EvaluateKids2 -- same as Evaluate Kids, but include VTK octree generation
336 //--------------------------------------------------------------------------------------
337 void Coupling::EvaluateKids2( const Vector3r& size, const int& level, const Vector3r& cpos,
338 const Complex& parentVal, vtkHyperOctree* oct, vtkHyperOctreeCursor* curse) {
339
340 std::cout << "\r" << (int)(1e2*VOLSUM/(Size[0]*Size[1]*Size[2])) << "\t" << nleaves;
341 std::cout.flush();
342
343 // Next level step, interested in one level below
344 // bitshift requires one extra, faster than, and equivalent to std::pow(2, level+1)
345 Vector3r step = size.array() / (Real)(1 << (level+1) );
346 Real vol = (step(0)*step(1)*step(2)); // volume of each child
347
348 Vector3r pos = cpos - step/2.;
349 Eigen::Matrix<Real, 8, 3> posadd = (Eigen::Matrix<Real, 8, 3>() <<
350 0, 0, 0,
351 step[0], 0, 0,
352 0, step[1], 0,
353 step[0], step[1], 0,
354 0, 0, step[2],
355 step[0], 0, step[2],
356 0, step[1], step[2],
357 step[0], step[1], step[2] ).finished();
358
359 VectorXcr kvals(8); // individual kernel vals
360 cpoints->ClearFields();
361 for (int ichild=0; ichild<8; ++ichild) {
362 Vector3r cp = pos; // Eigen complains about combining these
363 cp += posadd.row(ichild);
364 cpoints->SetLocation( ichild, cp );
365 }
366
367 Eigen::Matrix<Complex, 3, 8> Ht = Eigen::Matrix<Complex, 3, 8>::Zero();
368 Eigen::Matrix<Complex, 3, 8> Hr = Eigen::Matrix<Complex, 3, 8>::Zero();
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 oct->SubdivideLeaf(curse);
398 for (int ichild=0; ichild<8; ++ichild) {
399 curse->ToChild(ichild);
400 Vector3r cp = pos; // Eigen complains about combining these
401 cp += posadd.row(ichild);
402 /* Test for position via alternative means */
403 /*
404 Real p[3];
405 GetPosition(curse, p);
406 if ( (Vector3r(p) - cp).norm() > 1e-8 ) {
407 std::cout << "ERROR @ nleaves" << nleaves << "\n" << cp[0] << "\t" << p[0] << "\t" << cp[1] << "\t" << p[1]
408 << "\t" << cp[2] << "\t" << p[2] << "\t" << vol<< std::endl;
409 throw std::runtime_error("doom");
410 }
411 */
412 /* End of position test */
413 EvaluateKids2( size, level+1, cp, kvals.row(ichild), oct, curse );
414 curse->ToParent();
415 }
416 return; // not a leaf
417 }
418 LeafDict[curse->GetLeafId()] = ksum/(8.*vol);
419 LeafDictIdx[curse->GetLeafId()] = nleaves;
420 Kern.row(ilay) += ksum;
421 VOLSUM += 8*vol;
422 nleaves += 1;
423 return; // is a leaf
424 }
425
426 //--------------------------------------------------------------------------------------
427 // Class: Coupling
428 // Method: GetPosition
429 //--------------------------------------------------------------------------------------
430 void Coupling::GetPosition( vtkHyperOctreeCursor* Cursor, Real* p ) {
431 Real ratio=1.0/(1<<(Cursor->GetCurrentLevel()));
432 //step = ((Size).array() / std::pow(2.,Cursor->GetCurrentLevel()));
433 p[0]=(Cursor->GetIndex(0)+.5)*ratio*this->Size[0]+this->Origin[0] ;//+ .5*step[0];
434 p[1]=(Cursor->GetIndex(1)+.5)*ratio*this->Size[1]+this->Origin[1] ;//+ .5*step[1];
435 p[2]=(Cursor->GetIndex(2)+.5)*ratio*this->Size[2]+this->Origin[2] ;//+ .5*step[2];
436 }
437
438 #endif
439
440 } // ---- end of namespace Lemma ----
441
442 /* vim: set tabstop=4 expandtab */
443 /* vim: set filetype=cpp */
444