-
Notifications
You must be signed in to change notification settings - Fork 3
/
model.cpp
1228 lines (1052 loc) · 46.8 KB
/
model.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// model.cpp
//
#include <cstdlib> /* rand(), srand(), RAND_MAX */
#include <ctime> /* time() */
#include <iomanip>
#include "model.hpp"
#include "phonons.hpp"
#include "dataout.hpp"
// In this file:
// CLASS IMPLEMENTATIONS FOR:
//
// o Class ModelParams
// o Class Model
//
// Search on "&&&&" to jump between class implementations in this
// file.
//
//////////////////////////////////////////////////////////////////////////
// &&&& ****
// **** CLASS: ModelParams ****
// **** ****
//
//////
// METHOD: ModelParams :: AddSeismometerByWavelength()
//
void ModelParams::AddSeismometerByWavelength(
EarthCoords::Generic location,
axes_scheme_e orientation,
Real radius_wl) {
SeisRequest request;
request.Location = location;
request.Orientation = orientation;
request.GatherRadiusInner[RAY_P] = 0;
request.GatherRadiusInner[RAY_S] = 0;
request.GatherRadiusOuter[RAY_P] = radius_wl;
request.GatherRadiusOuter[RAY_S] = radius_wl;
request.RadiiUnitsAreWavelengths = true; // Signals that radii must
// be converted before
// constructing the
// Seismometer object
mSReqList.push_back(request); // Add to the request list
}
//////
// METHOD: ModelParams :: AddSeismometerFixedRadius()
//
void ModelParams::AddSeismometerFixedRadius(
EarthCoords::Generic location,
axes_scheme_e orientation,
Real radius ) {
SeisRequest request;
request.Location = location;
request.Orientation = orientation;
request.GatherRadiusInner[RAY_P] = 0;
request.GatherRadiusInner[RAY_S] = 0;
request.GatherRadiusOuter[RAY_P] = radius;
request.GatherRadiusOuter[RAY_S] = radius;
request.RadiiUnitsAreWavelengths = false;
mSReqList.push_back(request); // Add to the request list
}
//////
// METHOD: ModelParams :: AddSeismometerRing()
//
void ModelParams::AddSeismometerRing(
EarthCoords::Generic location,
axes_scheme_e orientation,
Real inner_radius,
Real outer_radius ) {
SeisRequest request;
request.Location = location;
request.Orientation = orientation;
request.GatherRadiusInner[RAY_P] = inner_radius;
request.GatherRadiusInner[RAY_S] = inner_radius;
request.GatherRadiusOuter[RAY_P] = outer_radius;
request.GatherRadiusOuter[RAY_S] = outer_radius;
request.RadiiUnitsAreWavelengths = false;
mSReqList.push_back(request); // Add to the request list
}
//////
// METHOD: ModelParams :: Output()
//
// Writes model parameters to stdout.
//
void ModelParams::Output() const
{
using std::cout;
using std::endl;
using std::setw;
ModelParams def;
std::cout
<< "TOA_Degree: " << TOA_Degree << setw(30)
<< ((TOA_Degree==def.TOA_Degree) ? "(default)" : "") << endl
<< "Frequency: " << Frequency << setw(31)
<< ((Frequency == def.Frequency) ? "(default)" : "") << endl
<< "Number of Phonons: " << NumPhonons << setw(22)
<< ((NumPhonons == def.NumPhonons) ? "(default)" : "") << endl
<< "Time to Live: " << PhononTTL << setw(27)
<< ((PhononTTL == def.PhononTTL) ? "(default)" : "") << endl
<< "Time Bin Size: " << GetBinSize() << setw(27)
<< ((GetBinSize() == def.GetBinSize()) ? "(default)" : "")
<< endl
<< "Cylinder Range: " << CylinderRange << setw(24)
<< ((CylinderRange == def.CylinderRange) ? "(default)" : "")
<< endl;
std::cout << "Args provided for \"compiled\" model are: ";
if (CompiledArgs.size() > 0) {
cout << "(nu,eps,a,k,Q)" << endl << " ";
for (Index i = 0; i < CompiledArgs.size(); i++) {
std::cout << CompiledArgs[i] << " ";
}
cout << endl;
} else {
cout << "(none)" << endl;
}
std::cout << mSReqList.size() << " Seismometers Requested." << endl;
std::cout
<< "Event Location: (" << EventSourceLoc.x1() << ", "
<< EventSourceLoc.x2() << ", " << EventSourceLoc.x3() << ") "
<< "(User Coord System)" << endl
<< "Event Moment Tensor:" << endl;
EventSourceMT.OutputContents();
}
//////
// METHOD: ModelParams :: OutputOctaveText()
//
// Writes the set of model parameters to an output stream (perhaps a
// file stream set up by the caller) in a text-based format that can
// be read by GNU/Octave and easily parsed in an octave scripts.
//
void ModelParams::OutputOctaveText(std::ostream * out) const {
*out << "# Generated by Radiative3D for input into GNU/Octave"
<< std::endl << std::endl;
*out << "# name: NumPhonons \n"
<< "# type: scalar \n"
<< NumPhonons << " \n" << " \n";
*out << "# name: TOA_Degree \n"
<< "# type: scalar \n"
<< TOA_Degree << " \n" << " \n";
*out << "# name: Frequency \n"
<< "# type: scalar \n"
<< Frequency << " \n" << " \n";
*out << "# name: PhononTTL \n"
<< "# type: scalar \n"
<< PhononTTL << " \n" << " \n";
*out << "# name: EventSourceLocUCS \n" // UCS means User-choice Coordinate
<< "# type: matrix \n" // System, as opposed to internal
<< "# rows: 1 \n" // XYZ model space coordinate
<< "# columns: 3 \n" // system.
<< EventSourceLoc.x1() << " "
<< EventSourceLoc.x2() << " "
<< EventSourceLoc.x3() << "\n" << "\n";
*out << "# name: EventSourceMT \n"
<< "# type: matrix \n"
<< "# rows: 3 \n"
<< "# columns: 3 \n"
<< EventSourceMT.xx() << " "
<< EventSourceMT.xy() << " "
<< EventSourceMT.xz() << " \n"
<< EventSourceMT.yx() << " "
<< EventSourceMT.yy() << " "
<< EventSourceMT.yz() << " \n"
<< EventSourceMT.zx() << " "
<< EventSourceMT.zy() << " "
<< EventSourceMT.zz() << " \n"
<< "\n";
*out << "# name: GridSource \n"
<< "# type: scalar \n"
<< GridSource << " \n" << " \n";
*out << "# name: CylinderRange \n"
<< "# type: scalar \n"
<< CylinderRange << " \n" << " \n";
*out << "# name: CompiledArgs \n"
<< "# type: matrix \n"
<< "# rows: 1 \n"
<< "# columns: " << CompiledArgs.size() << " \n";
for (Index i=0; i<CompiledArgs.size(); i++) {
*out << CompiledArgs[i] << " ";
}
*out << "\n\n";
}
//////////////////////////////////////////////////////////////////////////
// &&&& ****
// **** CLASS: Model ****
// **** ****
//
//////
// CONSTRUCTOR: ::: Model() :::
//
// Initializes our physical universe and computational environment,
// and builds our Earth model based on user wishes as communicated to
// us via a ModelParams object.
//
// Generally, only ONE model object should be instantiated at any
// given time. Among the responsibilities of this constructor is to
// initialize certain global and class-static parameters, and
// instantiating multiple Model objects could result in them
// stepping on each others' toes.
//
//
Model::Model(const ModelParams & par) {
std::cout << "@@ __BEGIN_MODEL_INITIALIZATION__" << std::endl;
// ::::::
// :: Initialize members of the Model object:
// :
mNumPhonons = par.NumPhonons;
// ::::::
// :: Initialize global standard-lib objects:
// :
srand(time(NULL)); // Seed the random number generator.
// (TODO: Look into better random number
// implementations than the standard library
// RNG.)
// ::::::
// :: Create objects that will be used/needed by other classes:
// :
//
// Create/Initialize the TakeOffAngle array,
// which is a discrete set of precomputed
// directions used by event-source and
// scatterer PhononSource objects:
//
std::cout << "@@ __DISCRETIZING_TOA__" << std::endl;
mpTOA_Array = new S2::TesselSphere(S2::TESS_ICO,par.TOA_Degree);
std::cout << "|\n";
std::cout
<< "|" << std::setw(8) << mpTOA_Array->size()
<< " Take-off angles initialized for event and scattering sources.\n";
std::cout
<< "| (TesselSphere of degree " << par.TOA_Degree << ".)\n";
std::cout << "|\n";
// ::::::
// :: Inform other classes of needed parameters:
// :
//
// The PhononSource class needs to know about the take-off
// angles array:
//
PhononSource::Set_TOA_Array(mpTOA_Array);
//
// The Phonon class needs to know the maximum sim-time to
// compute per phonon before abandoning it as irrelevant:
//
Phonon::SetTimeToLive(par.PhononTTL);
//
// The ScatterParams class and the MediumCell class need
// to know the chosen phonon frequency (ScatterParams uses
// it to automatically compute 'el' and 'gama0' from
// velocity parameters, and MediumCell uses it to compute
// wavelengths from velocities):
//
ScatterParams::SetFrequencyHertz(par.Frequency);
MediumCell ::SetFrequencyHertz(par.Frequency);
//
// The Seismometer class needs to know the desired width
// of time bins and the desired recording duration:
//
Seismometer::SetTimeBinParameters( //
par.GetBinSize(), // Bin size
par.PhononTTL // Recording time
); //
// ::::::
// :: Build the Earth model:
// :
//
// We start by building the Grid - the lattice of nodes
// that establish the geometry and physical properties of
// the earth model.
//
std::cout << "@@ __CONSTRUCTING_GRID__" << std::endl;
std::cout << "|\n";
switch ( par.GridSource ) {
case ModelParams::GRID_UNSPEC: // Grid source unspecified
std::cout << "| Grid source unspecified. " // TEMP
<< "Assuming user-compiled with index 0.\n"; // TEMP
mGrid.ConstructGridManual(0, par.CompiledArgs); // TEMP
break;
// TODO: Assume a "Built-In" test model instead.
case ModelParams::GRID_COMPILED: // Grid hard-coded by user
std::cout << "| Grid source: User-coded grid; Selection ID = "
<< par.CompiledSelector << " with "
<< par.CompiledArgs.size() << " args.\n";
mGrid.ConstructGridManual( // (defined in user.cpp)
par.CompiledSelector, par.CompiledArgs);
break;
case ModelParams::GRID_FROMFILE: // Grid comes from file
// TEMP: Fall-through to default (not implemented yet)
// ::::::
// :: Grid-Type: Unrecognized
// :
default:
TextStream NoGridHandlerMsg;
NoGridHandlerMsg << "No handler exists for Grid Source '"
<< par.GridSource << "'.";
throw(Runtime(NoGridHandlerMsg.str()));
break;
} // END switch (par.GridSource)
///
std::cout << "|\n"; // Tell user what we did:
std::cout << "| " << std::setw(6) << mGrid.N()
<< " Grid nodes initialized. Arrangement: "
<< mGrid.Ni() << " x " << mGrid.Nj() << " x "
<< mGrid.Nk() << "\n|\n";
if (!ECS.IsEarthFlattening() && !ECS.CurvedCoords()) std::cout
<< "| Coordinate mapping is rectilinear.\n";
if (ECS.IsEarthFlattening()) std::cout
<< "| An Earth-flattening transformation was used with"
<< " Earth Radius = " << ECS.GetEarthRadius() << "\n";
if (ECS.CurvedCoords()) std::cout
<< "| Level planes were curved with Earth Radius = "
<< ECS.GetEarthRadius() << "\n";
std::cout << "|\n";
//
// Next we "fill in" the grid with "cells" that fill the
// volume of space, modelling elastic media according to
// the values known at the vertices. We have a few
// different cell-types to choose from, all of which are
// derivatives of the MediumCell virtual base class.
//
std::cout << "@@ __STAGING_EARTH_MODEL_ONTO_GRID__" << std::endl;
std::cout << "|\n";
Text WaitNotice = (par.TOA_Degree > 7) // With lots of TOAs alloc'ing
? "| (This may take a while.)\n" : ""; // the scatterers takes a long
// time...
switch ( mGrid.GetModelType() ) {
case Grid::MOD_CYLINDER: // Layered aka Stacked-Cylinder Model:
//
std::cout << "| Building tilted-interface LAYERED model"
<< " from plumbline grid... \n"
<< WaitNotice << std::flush;
RCUCylinder::SetRange(par.CylinderRange);
BuildCellArray_Cylinder();
break;
case Grid::MOD_TETRAWCG: // Tetrahedral model from WCG Grid:
//
std::cout << "| Building TETRA model from Warped Cartesian Grid... \n"
<< WaitNotice << std::flush;
BuildCellArray_WCGTetra();
break;
case Grid::MOD_SPHERESHELL: // Spherical shell model from dropline Grid:
//
std::cout << "| Building SPHERICAL SHELL model... \n"
<< WaitNotice << std::flush;
BuildCellArray_SphericalShells();
break;
default:
TextStream NoModelHandlerMsg;
NoModelHandlerMsg << "No handler exists for Model type '"
<< mGrid.GetModelType() << "'.";
throw(Runtime(NoModelHandlerMsg.str()));
break;
}
std::cout << "|\n"; // Output model cell statististics
std::cout
<< "|" << std::setw(8) << mCellArray.size()
<< " Model cells constructed.\n";
std::cout
<< "|" << std::setw(8) << Scatterer::GetCollectionSize()
<< " Scatterer objects allocated among cells.\n";
std::cout << "|\n";
//
// Finally, we create a PhononSource object to represent
// the seismic event. The object will be of type
// ShearDislocation and will be defined by the moment-
// tensor argument provided in the parameters object:
//
std::cout << "@@ __INITIALIZING_EVENT_SOURCE__" << std::endl;
Tensor::Tensor eventMT = par.EventSourceMT;
R3::XYZ eventLoc = ECS.Convert(par.EventSourceLoc);
eventMT.Transform(ECS.GetXYZToLocalNEDRotation(eventLoc));
// Assume the moment tensor that has been provided
// uses the xyz <--> NED (North, East, Down)
// coordinate system, which might not be in agreement
// with the local orientation of the Earth model at
// the event source location. Rotate the moment
// tensor so the the user gets the expected result.
mpEventSource = new ShearDislocation(eventMT, eventLoc);
mpEventSource->Link(FindCellContainingPoint(eventLoc));
// Create EventSource and inform it of
// the MediumCell in which it resides.
Seismometer::SetEventParameters(eventLoc, par.EventSourceMT);
// Inform seismometer class of event parameters
// (used to output metadata with traces.)
// ::::::
// :: Create the Seismometers:
// :
std::cout << "@@ __INITIALIZING_SEISMOMETERS__" << std::endl;
for ( std::vector<ModelParams::SeisRequest>::const_iterator
sr = par.mSReqList.begin();
sr != par.mSReqList.end();
++sr) {
// FOR EACH SEISREQUEST "sr"; DO:
//
const R3::XYZ Location = FindSurface(ECS.Convert(sr->Location));
Real rad_in [RAY_NUMBASICTYPES]; // Inner-Radii
Real rad_out[RAY_NUMBASICTYPES]; // Outer-Radii
// :: Convert radii units if needed:
if (sr->RadiiUnitsAreWavelengths == false)
// Then no conversion is needed:
{
rad_in[RAY_P] = sr->GatherRadiusInner[RAY_P];
rad_in[RAY_S] = sr->GatherRadiusInner[RAY_S];
rad_out[RAY_P] = sr->GatherRadiusOuter[RAY_P];
rad_out[RAY_S] = sr->GatherRadiusOuter[RAY_S];
} else {
// Then conversion from wavelengths
// to model units is needed:
Real wavlen[RAY_NUMBASICTYPES]; // Wavelengths
MediumCell * pCC; // ContainerCell
pCC = FindCellContainingPoint(Location);
wavlen[RAY_P] = pCC->GetWavelengthAtPoint(Location, RAY_P);
wavlen[RAY_S] = pCC->GetWavelengthAtPoint(Location, RAY_S);
rad_in [RAY_P] = sr->GatherRadiusInner[RAY_P] * wavlen[RAY_P];
rad_in [RAY_S] = sr->GatherRadiusInner[RAY_S] * wavlen[RAY_S];
rad_out[RAY_P] = sr->GatherRadiusOuter[RAY_P] * wavlen[RAY_P];
rad_out[RAY_S] = sr->GatherRadiusOuter[RAY_S] * wavlen[RAY_S];
}
R3::XYZ SeisAxisX1 = ECS.GetEast(Location);
std::string AxesDesc = "ENZ";
if (sr->Orientation == ModelParams::AX_RTZ) {
SeisAxisX1 = ECS.GetRadial(eventLoc,Location);
AxesDesc = "RTZ";
}
// :: Ask the DataOut object to actually create the Seismometer:
dataout.AddSeismometer(Location, SeisAxisX1, rad_in, rad_out, AxesDesc);
//
// END (FOR EACH SEISREQUEST)
}//
std::cout << "@@ __MODEL_INITIALIZATION_COMPLETE__"
<< std::endl << std::flush;
}
Model::~Model() {
delete mpEventSource;
delete mpTOA_Array;
}
//////
// METHOD: FindCellContainingPoint()
//
// Walks the array of MediumCells and returns a pointer to the one
// that "best" contains the R3 point location provided in the
// arguments. "Best" here refers to the fact that a point
// intended to be "just on the surface" of the cell (e.g., the
// location of a seismometer on the Earth's surface) might, due to
// round-off error, calculate as being slightly outside the
// cell. This function will return that cell anyway. If no cell
// adequately contains the point, a null pointer is returned.
//
MediumCell * Model::FindCellContainingPoint(const R3::XYZ & loc) const {
Real best_mismatch; // Goal is to "minimize" the mismatch
MediumCell * best_mcellp; //
best_mismatch = 0.1; // Max tolerable mismatch (100 meters)
best_mcellp = 0; // Return a default null ptr if no tolerable
// cell is found
MediumCellPtrArray::const_iterator cellp_it = mCellArray.begin();
MediumCellPtrArray::const_iterator cellp_end = mCellArray.end();
for (; cellp_it != cellp_end; cellp_it++) {
Real mismatch = (*cellp_it)->IsPointInside(loc);
if (mismatch <= 0.) {
best_mismatch = mismatch; // Then we're definitively inside
best_mcellp = *cellp_it; // Return this cell
break; //
}
else if (mismatch < best_mismatch) {
best_mismatch = mismatch; // Then it's a better match than previous
best_mcellp = *cellp_it; // but keep searching
}
}
return best_mcellp;
}
//////
// METHOD: Model :: FindSurface()
//
// Given a location presumed to be near the surface, this method
// finds a corresponding point on the surface that is found by
// following a vertical (skyward) line that intersects both the
// point and the surface.
//
R3::XYZ Model::FindSurface(R3::XYZ loc) const {
if (mSurfaceFaces.size() > 1) { // TODO: write flexible version of this
// For multiple surface-face models (e.g. in a tetrahedral model),
// the process of surface-finding is more complex, ...and hasn't
// been written yet.
throw(std::invalid_argument(
"FindSurface() doesn't know how to handle multiple surface faces yet."));
return loc;
} else if (mSurfaceFaces.size() == 0) {
// If the model-building code didn't record any surface faces
// (possible for test models I guess) then we silently return the
// starting location. TODO: Ponder whether a warning should be
// issued for this.
return loc;
} // Else proceed:
CellFace & sface = *mSurfaceFaces[0];
R3::XYZ skyward = ECS.GetUp(loc);
Real dist = sface.LinearRayDistToExit(loc,skyward);
R3::XYZ new_loc = loc + skyward.ScaledBy(dist);
if (false) { // (set true for debug output)
std::cout << "Nudged surface location from ("
<< loc.x() << ", " << loc.y() << ", " << loc.z()
<< ") to ("
<< new_loc.x() << ", " << new_loc.y() << ", " << new_loc.z()
<< ") dist was: " << dist << " sface was at: " << &sface << "\n";
}
return new_loc;
}
//////
// METHOD: Model :: RunSimulation()
//
// Runs the simulation outer-loop.
//
void Model::RunSimulation() {
unsigned nph_byten = mNumPhonons / 10; // For Progress Indicator
unsigned nph_by100 = mNumPhonons / 100; //
bool finegrain = (mNumPhonons > 9999); // Threshold for hundredths indicator
if (nph_byten==0) nph_byten=1; // Prevent div/0
std::cout << "@@ __BEGINNING_SIMULATION__" << std::endl << std::flush;
for (int i=0; i < mNumPhonons; i++) {
Phonon P = mpEventSource->GenerateEventPhonon();
P.Propagate();
if ( finegrain && (i % nph_by100 == 0) && (i>0) ) { // 1% mark
std::cerr << "."; //
}
if ( i % nph_byten == 0 ) { // 10% mark
unsigned pct = (i / nph_byten) * 10; //
std::cerr << pct << "% of " << mNumPhonons
<< " have been cast.\n";
}
}
if (finegrain) {std::cerr << ".";}
std::cerr << "100% of " << mNumPhonons << " have been cast.\n";
std::cout << "@@ __SIMULATION_COMPLETE__" << std::endl;
dataout.OutputPostSimSummary();
}
//////
// METHOD: Model :: BuildCellArray_Cylinder()
//
// The BuildCellArray family of helper functions assist in the
// construction of the earth model by building the array of
// MediumCell objects that make up the material/volumetric part of
// the model. This one builds the cell array out of Cylinder
// MediumCell types. The information needed to build the array
// comes from the Grid object (mGrid), which must be complete before
// this function is called.
//
void Model::BuildCellArray_Cylinder() {
// ::::::
// :: Reserve space in the cell array:
// :
// (Can speed allocation for very large models. Of
// course, we're only reserving room for pointers
// here, but still.)
//
unsigned int array_num_cells;
array_num_cells = mGrid.Nk() - 1; // Pretty simple calculation for
// the cylinder model, which is
// just a vertical stack.
mCellArray.reserve(array_num_cells);
// ::::::
// :: Walk the grid from top-to-bottom (along k-index) and create
// :: MediumCells along the way:
// :
for (unsigned k = 0; k < array_num_cells; k++) {
// GridData: Cell is bounded on top by the underside of grid
// layer k, and on bottom by the above-side grid
// layer k+1
const GridData GD_Top = mGrid.Node(0,0,k).Data(GridNode::GN_BELOW);
const GridData GD_Bot = mGrid.Node(0,0,k+1).Data(GridNode::GN_ABOVE);
mCellArray.push_back( // Construct the cell
new RCUCylinder(mGrid.Node(0,0,k).Loc(),
mGrid.Node(1,0,k).Loc(),
mGrid.Node(2,0,k).Loc(),
mGrid.Node(0,0,k+1).Loc(),
mGrid.Node(1,0,k+1).Loc(),
mGrid.Node(2,0,k+1).Loc(),
GD_Top.getV(), GD_Bot.getV(),
GD_Top.Rho(), GD_Bot.Rho(),
GD_Top.getQ(), GD_Bot.getQ()));
ScatterParams scatpar(GD_Top.getV(), GD_Top.getHS());
Scatterer * pscat = Scatterer::GetScattererMatchingParams(scatpar);
mCellArray[k]->Link(pscat);
}// end for(k)
//
// ::::::
// :: Walk the Cylinder-stack from top to bottom and link the
// :: interfaces between them along the way:
// :
for (unsigned k = 1; k < array_num_cells; k++) {
// CellFace References:
// upper: the BOTTOM face of the cell ABOVE
// lower: the TOP face of the cell BELOW
CellFace & upper = mCellArray[k-1]->Face(CellFace::F_BOTTOM);
CellFace & lower = mCellArray[k] ->Face(CellFace::F_TOP);
bool dis = mGrid.Node(0,0,k).IsDiscontinuous();
upper.LinkTo(lower, dis); // Lets the faces know about each
// other, and sets the appropriate
// connectivity flags
}
// ::::::
// :: Set collection and reflection status at the top boundary
// :: of the model:
// :
CellFace & ModelTop = mCellArray[0]->Face(CellFace::F_TOP);
//CellFace & ModelBot = mCellArray.back()->Face(CellFace::F_BOTTOM);
ModelTop.SetCollect(true);
ModelTop.SetReflect(true);
mSurfaceFaces.push_back(&ModelTop); // Add to surface face list,
// for surface-finding routine
// in FindSurface().
}
//////
// METHOD: Model :: BuildCellArray_SphericalShells()
//
// Builds out the cell array using SphereShell cell types.
//
void Model::BuildCellArray_SphericalShells() {
// Reserve space in the cell array:
unsigned int array_num_cells;
array_num_cells = mGrid.Nk() - 1;
mCellArray.reserve(array_num_cells);
// ::::::
// :: Walk the grid from top-to-bottom (along k-index) and create
// :: MediumCells along the way:
// :
for (unsigned k = 0; k < array_num_cells; k++) {
// GridData: Cell is bounded on top by the underside of grid
// layer k, and on bottom by the above-side grid
// layer k+1
const GridData GD_Top = mGrid.Node(0,0,k).Data(GridNode::GN_BELOW);
const GridData GD_Bot = mGrid.Node(0,0,k+1).Data(GridNode::GN_ABOVE);
mCellArray.push_back( // Construct the cell
new SphereShell(mGrid.Node(0,0,k).GetRawLoc().Radius(ECS),
mGrid.Node(0,0,k+1).GetRawLoc().Radius(ECS),
GD_Top, GD_Bot)
);
ScatterParams scatpar(GD_Top.getV(), GD_Top.getHS());
Scatterer * pscat = Scatterer::GetScattererMatchingParams(scatpar);
mCellArray[k]->Link(pscat);
}// end for(k)
//
// ::::::
// :: Walk the Cylinder-stack from top to bottom and link the
// :: interfaces between them along the way:
// :
for (unsigned k = 1; k < array_num_cells; k++) {
// CellFace References:
// upper: the BOTTOM face of the cell ABOVE
// lower: the TOP face of the cell BELOW
CellFace & upper = mCellArray[k-1]->Face(CellFace::F_BOTTOM);
CellFace & lower = mCellArray[k] ->Face(CellFace::F_TOP);
bool dis = mGrid.Node(0,0,k).IsDiscontinuous();
upper.LinkTo(lower, dis); // Lets the faces know about each
// other, and sets the appropriate
// connectivity flags
}
// ::::::
// :: Set collection and reflection status at the top boundary
// :: of the model:
// :
CellFace & ModelTop = mCellArray[0]->Face(CellFace::F_TOP);
//CellFace & ModelBot = mCellArray.back()->Face(CellFace::F_BOTTOM);
ModelTop.SetCollect(true);
ModelTop.SetReflect(true);
mSurfaceFaces.push_back(&ModelTop); // Add to surface face list,
// for surface-finding routine
// in FindSurface().
}
//////
// METHOD: Model :: BuildCellArray_WCGTetra()
//
// Builds out the cell array, using tetrahedral cell, assuming the
// grid structure is "Warped Cartesian Grid." This means that the
// grid codes for cuboidal subunits in which properties are
// specified on the eight corners, and the subunits are arranged in
// a regular array of rows, columns, and stacks. Each subunit will
// comprise five individual tetra cells. Note that it is up to the
// user how the row indx (i), column index (j), and stack index (k)
// map to real-life spatial coordinates (XYZ), but the assumption
// made for model building is that k increases with increasing
// depth. This is only important when the user wants discontinuous
// interfaces, where there is a "sudden" change in property values
// "above" (in direction of smaller k) and "below" (in direction of
// larger k) a particular interface.
//
// RESPONSIBILITIES: Builds array of tetra cells. Links them
// together within a cuboidal subunit (actually BuildBasicPattern
// handles this), and accross cuboidal subunits (via calls to
// LinkBlocks). And ensures that surface interface is marked
// "collecting" and "reflecting".
//
// GRID INTERPRETATION: Nodes be like: (ijk indices)
//
// 000-------010-------020------- .
// |\ \ \ .
// | \ \ \ o----> (+j) Northwards (typical) .
// | \ \ \ |\ .
// | 100-------110-------120------- | \ .
// | |\ |\ |\ | (+i) Eastwards (typical) .
// 001 | \ | \ | \ v .
// |\ | \ | \ | \ (+k) Downwards (typical) .
// | \ | | | .
// | \| | | .
// | 101-------111-------121------- .
// | |\ |\ |\ .
// | \ | \ | \ .
//
//
// BUILDOUT ORDER: Subunits are built in plunging order first (along
// k index), then what we call "row-building" order (along j index),
// and finally "column-building" order (along i index). Adhering to
// this order means that the five tetras comprising a particular ijk
// subunit will have base index into the cell array of:
//
// base = ((i*(nK*nJ) + j*nK + k) * 5)
//
// this formula is codified in helper function WCGBaseIndexFromIJK(),
// which should be used (rather than coding the formula directly)
// whenever a base index is needed in order to reduce code duplication
// (and resultant opportunities for coding error).
//
// CELL ORDER WITHIN SUBUNIT: There are five tetra cells per
// subunit. The internal tetra (touches no exterior faces of
// subunit) of a particular subunit is always at offset 0 from the
// base index. The tetra at offset 1 is the one covering the
// southwest (small i, small j) edge of the subunit. Offset 1
// covers the northwest corner. Offset 2 covers the northeast
// corner. And offset 4 covers the southeast corner.
//
// MIRROR PARITY: In order for the cell faces across subunits to
// join up correctly, the basic pattern must invert via a reflection
// as one moves from one subunit to an adjacent one. The basic
// pattern and its reflection is depicted below, and is described in
// further detail in the commentary to WCGBuildBasicPattern() and
// WCGLinkBlocksForward().
//
// 000-------010-------020 .
// |\ * (2) \ (1) * \ .
// |.\ * \ * \ o----> (+j) Northwards (typical) .
// | \ (4) *\ * (3) \ |\ .
// | 100-------110-------120 | \ .
// | . | *|* | | (+i) Eastwards (typical) .
// 001 | (4) * | * (3) | v .
// \ | * | * | (+k) Downwards (typical) .
// \'| * (3) | (4) * | .
// \|* | *| Tetra offsets labeled as (n) .
// 101-------111-------121 .
// Natural Reflected .
//
//
void Model::BuildCellArray_WCGTetra() {
int nI = mGrid.Ni()-1; // Number of rows (of pattern blocks)
int nJ = mGrid.Nj()-1; // Number of columns
int nK = mGrid.Nk()-1; // Number of stacks
mCellArray.clear();
mCellArray.reserve(nI*nJ*nK*5);
// Vector of grid nodes used to construct 'basic pattern' blocks.
// Each block is constructed using 8 grid nodes.
// (Adjacent blocks share grid nodes.)
std::vector<const GridNode *> nodes;
nodes.reserve(8); // Room for an octuple of nodes.
// Build all blocks:
for (int i = 0; i < nI; i++) {
for (int j = 0; j < nJ; j++) {
for (int k = 0; k < nK; k++) {
// Build out node octuple
nodes.clear(); // (only one octuple needed per basic pattern)
nodes.push_back(&mGrid.RelNode(i,j,k,0,0,0));
nodes.push_back(&mGrid.RelNode(i,j,k,0,0,1));
nodes.push_back(&mGrid.RelNode(i,j,k,0,1,0));
nodes.push_back(&mGrid.RelNode(i,j,k,0,1,1));
nodes.push_back(&mGrid.RelNode(i,j,k,1,0,0));
nodes.push_back(&mGrid.RelNode(i,j,k,1,0,1));
nodes.push_back(&mGrid.RelNode(i,j,k,1,1,0));
nodes.push_back(&mGrid.RelNode(i,j,k,1,1,1));
// Build block
bool mirror = (((i+j+k)%2)==1); // Reflection tracking
bool issurface = (k==0); // Surface is top of stacks
WCGBuildBasicPattern(nodes,mirror,issurface);
// Link Backwards:
Index Block = WCGBaseIndexFromIJK(i,j,k,nI,nJ,nK); // This block
bool adjmirror = !mirror; // If 'Block' is reflected, then
// adjacent blocks are not.
if (i>0) { // Then there's a block behind us:
Index Previous = WCGBaseIndexFromIJK(i-1,j,k,nI,nJ,nK);
WCGLinkBlocksForward(Previous,Block,CellFace::FACE_C,adjmirror);
}
if (j>0) { // Then there's a block to the left of us:
Index Previous = WCGBaseIndexFromIJK(i,j-1,k,nI,nJ,nK);
WCGLinkBlocksForward(Previous,Block,CellFace::FACE_B,adjmirror);
}
if (k>0) { // Then there's a block above us:
Index Previous = WCGBaseIndexFromIJK(i,j,k-1,nI,nJ,nK);
WCGLinkBlocksForward(Previous,Block,CellFace::FACE_D,adjmirror);
}
}
}
}// End triple for Loop
// All blocks now built and linked.
}// End BuildCellArray_WCGTetra()
//
//////
// METHOD Model :: WCGBaseIndexFromIJK() [static]
//
// Helper function that computes the base index (index of first of
// five Tetra forming a block) of the i,j,k block in the cell array
// that results from Warped Cartesian (WCG) interpretation of the
// Grid array. The formula that computes this base index is given
// in the commentary for BuildCellArray_WCGTetra() and conforms to
// the build-out order of blocks during model construction. It is
// codified here to reduce code duplication.
//
Index Model::WCGBaseIndexFromIJK(Index i, Index j, Index k,
Count nI, Count nJ, Count nK) {
return ((i*(nK*nJ) + j*nK + k) * 5);
}
//////
// METHOD: Model :: WCGBuildBasicPattern()
//
// Constructs a "block" of five Tetra cells spanning a set of eight
// Grid Nodes representing a 2x2x2 subunit of a Warped Cartesian
// Grid (WCG) structured Grid collection. Pushes the five new Tetra
// cells onto the back end of mCellArray. Responsible for internal
// face linkages among the cells, (but NOT responsible for external
// face linkages between this and adjacent blocks). Responsible for
// requesting and linking Scatterers for each cell. Responsible for
// marking top surface and bottom surface faces as "discontinous" if
// indicated by the appropriate grid nodes. Responsible, if
// isSurface==true, for flagging the top surface faces as
// "collecting" and "reflecting.
//
// Note that the basic pattern has two variants, which we call
// "natural" and "reflected." Which pattern to build is a function
// of the mirror parameter (true --> build reflected pattern). An
// alternating pattern of reflected/non-reflected blocks means that
// they can be joined continuously across matched-up faces.
//
// Grid node references are passed in an eight element array. The
// pack order is illustrated by the square-bracketted numbers [n] in
// the figure below. The Tetra are constructed in a sequential
// order starting with the central (not visible) Tetra, and
// continuing in a clockwise fashion from there. The Tetra are
// indexed in the figure with parenthetical numbers (n). Each Tetra
// has four faces labeled A,B,C,D, and are oriented such that face A
// joins the internal Tetra, face B is flush with the j-index axis,
// face C with the i-index axis, and face D with the k-index axis.
// Face labels are also illustrated in the figure below.
//
//
// [0]-------------[2] [0]-------------[2]
// |\ * D \ |\ D * \ o------> (+j)
// |*\ * (2) \ | \ (1) * \ |\.
// |* \ * \ | \ * \ | \.
// | * \ (4) * \ |(1)\ * (3) \ | (+i)
// | * \ D * \ | B \ * D \ V
// |B * [4]-------------[6] | .[4]-------------[6] (+k)
// (1)*B | * | | .' | * |
// [1] * | (4) * | [1] B | * C |
// \ * | C * | \ (4)| * (3) |
// \ *| * | \ | * |
// \ *| * C | \ | (4) * |
// \*| * (3) | \ | C * |