forked from BLAST-WarpX/warpx
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBTDiagnostics.cpp
1504 lines (1374 loc) · 78.6 KB
/
BTDiagnostics.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
/* Copyright 2021 Revathi Jambunathan
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#include "BTDiagnostics.H"
#include "BTD_Plotfile_Header_Impl.H"
#include "ComputeDiagFunctors/BackTransformFunctor.H"
#include "ComputeDiagFunctors/CellCenterFunctor.H"
#include "ComputeDiagFunctors/ComputeDiagFunctor.H"
#include "ComputeDiagFunctors/RhoFunctor.H"
#include "Diagnostics/Diagnostics.H"
#include "Diagnostics/FlushFormats/FlushFormat.H"
#include "ComputeDiagFunctors/BackTransformParticleFunctor.H"
#include "Utils/Algorithms/IsIn.H"
#include "Utils/Parser/ParserUtils.H"
#include "Utils/TextMsg.H"
#include "Utils/WarpXConst.H"
#include "WarpX.H"
#include <ablastr/coarsen/sample.H>
#include <ablastr/utils/Communication.H>
#include <ablastr/utils/SignalHandling.H>
#include <ablastr/warn_manager/WarnManager.H>
#include <AMReX.H>
#include <AMReX_Algorithm.H>
#include <AMReX_BLassert.H>
#include <AMReX_BoxArray.H>
#include <AMReX_Config.H>
#include <AMReX_CoordSys.H>
#include <AMReX_DistributionMapping.H>
#include <AMReX_FileSystem.H>
#include <AMReX_ParallelContext.H>
#include <AMReX_ParallelDescriptor.H>
#include <AMReX_ParmParse.H>
#include <AMReX_Utility.H>
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <memory>
#include <sstream>
#include <vector>
using namespace amrex::literals;
namespace
{
constexpr int permission_flag_rwxrxrx = 0755;
}
BTDiagnostics::BTDiagnostics (int i, std::string name)
: Diagnostics(i, name)
{
ReadParameters();
}
void BTDiagnostics::DerivedInitData ()
{
auto & warpx = WarpX::GetInstance();
m_gamma_boost = WarpX::gamma_boost;
m_beta_boost = std::sqrt( 1._rt - 1._rt/( m_gamma_boost * m_gamma_boost) );
m_moving_window_dir = WarpX::moving_window_dir;
// Currently, for BTD, all the data is averaged+coarsened to coarsest level
// and then sliced+back-transformed+filled_to_buffer.
// The number of levels to be output is nlev_output.
nlev_output = 1;
m_t_lab.resize(m_num_buffers);
m_snapshot_domain_lab.resize(m_num_buffers);
m_buffer_domain_lab.resize(m_num_buffers);
m_snapshot_box.resize(m_num_buffers);
m_buffer_box.resize(m_num_buffers);
m_current_z_lab.resize(m_num_buffers);
m_current_z_boost.resize(m_num_buffers);
m_old_z_boost.resize(m_num_buffers);
m_buffer_counter.resize(m_num_buffers);
m_snapshot_ncells_lab.resize(m_num_buffers);
m_cell_centered_data.resize(nmax_lev);
m_cell_center_functors.resize(nmax_lev);
m_max_buffer_multifabs.resize(m_num_buffers);
m_buffer_flush_counter.resize(m_num_buffers);
m_geom_snapshot.resize( m_num_buffers );
m_snapshot_full.resize( m_num_buffers );
m_lastValidZSlice.resize( m_num_buffers );
m_buffer_k_index_hi.resize(m_num_buffers);
m_first_flush_after_restart.resize(m_num_buffers);
m_snapshot_geometry_defined.resize(m_num_buffers);
m_field_buffer_multifab_defined.resize(m_num_buffers);
for (int i = 0; i < m_num_buffers; ++i) {
m_geom_snapshot[i].resize(nmax_lev);
m_snapshot_full[i] = 0;
m_lastValidZSlice[i] = 0;
m_buffer_flush_counter[i] = 0;
m_first_flush_after_restart[i] = 1;
m_snapshot_geometry_defined[i] = 0;
m_field_buffer_multifab_defined[i] = 0;
}
for (int lev = 0; lev < nmax_lev; ++lev) {
// Define cell-centered multifab over the whole domain with
// user-defined crse_ratio for nlevels
DefineCellCenteredMultiFab(lev);
}
/* Allocate vector of particle buffer vectors for each snapshot */
MultiParticleContainer& mpc = warpx.GetPartContainer();
// If not specified, and write species is not 0, dump all species
const amrex::ParmParse pp_diag_name(m_diag_name);
int write_species = 1;
pp_diag_name.query("write_species", write_species);
if ((m_output_species_names.empty()) && (write_species == 1)) {
m_output_species_names = mpc.GetSpeciesNames();
}
m_do_back_transformed_particles =
((!m_output_species_names.empty()) && (write_species == 1));
// Turn on do_back_transformed_particles in the particle containers so that
// the tmp_particle_data is allocated and the data of the corresponding species is
// copied and stored in tmp_particle_data before particles are pushed.
if (m_do_back_transformed_particles) {
mpc.SetDoBackTransformedParticles(m_do_back_transformed_particles);
for (auto const& species : m_output_species_names){
mpc.SetDoBackTransformedParticles(species, m_do_back_transformed_particles);
}
}
m_particles_buffer.resize(m_num_buffers);
m_totalParticles_in_buffer.resize(m_num_buffers);
// check that simulation can fill all BTD snapshots
const int lev = 0;
const amrex::Real dt_boosted_frame = warpx.getdt(lev);
const int moving_dir = WarpX::moving_window_dir;
const amrex::Real Lz_lab = warpx.Geom(lev).ProbLength(moving_dir) / WarpX::gamma_boost / (1._rt+WarpX::beta_boost);
const int ref_ratio = 1;
const amrex::Real dz_snapshot_grid = dz_lab(dt_boosted_frame, ref_ratio);
// Need enough buffers so the snapshot length is longer than the lab frame length
// num_buffers * m_buffer_size * dz_snapshot_grid >= Lz
const int num_buffers = static_cast<int>(std::ceil(Lz_lab / m_buffer_size / dz_snapshot_grid));
const int final_snapshot_iteration = m_intervals.GetFinalIteration();
// the final snapshot starts filling when the
// right edge of the moving window intersects the final snapshot
// time of final snapshot : t_sn = t0 + i*dt_snapshot
// where t0 is the time of first BTD snapshot, t0 = zmax / c * beta / (1-beta)
//
// the right edge of the moving window at the time of the final snapshot
// has space time coordinates
// time t_intersect = t_sn, position z_intersect=zmax + c*t_sn
// the boosted time of this space time pair is
// t_intersect_boost = gamma * (t_intersect - beta * z_intersect_boost/c)
// = gamma * (t_sn * (1 - beta) - beta * zmax / c)
// = gamma * (zmax*beta/c + i*dt_snapshot*(1-beta) - beta*zmax/c)
// = gamma * i * dt_snapshot * (1-beta)
// = i * dt_snapshot / gamma / (1+beta)
//
// if j = final snapshot starting step, then we want to solve
// j dt_boosted_frame >= t_intersect_boost = i * dt_snapshot / gamma / (1+beta)
// j >= i / gamma / (1+beta) * dt_snapshot / dt_boosted_frame
const int final_snapshot_starting_step = static_cast<int>(std::ceil(final_snapshot_iteration / WarpX::gamma_boost / (1._rt+WarpX::beta_boost) * m_dt_snapshots_lab / dt_boosted_frame));
const int final_snapshot_fill_iteration = final_snapshot_starting_step + num_buffers * m_buffer_size - 1;
const amrex::Real final_snapshot_fill_time = final_snapshot_fill_iteration * dt_boosted_frame;
if (WarpX::compute_max_step_from_btd) {
if (final_snapshot_fill_iteration > warpx.maxStep()) {
warpx.updateMaxStep(final_snapshot_fill_iteration);
amrex::Print()<<"max_step insufficient to fill all BTD snapshots. Automatically increased to: "
<< final_snapshot_fill_iteration << std::endl;
}
if (final_snapshot_fill_time > warpx.stopTime()) {
warpx.updateStopTime(final_snapshot_fill_time);
amrex::Print()<<"stop_time insufficient to fill all BTD snapshots. Automatically increased to: "
<< final_snapshot_fill_time << std::endl;
}
if (warpx.maxStep() == std::numeric_limits<int>::max() && warpx.stopTime() == std::numeric_limits<amrex::Real>::max()) {
amrex::Print()<<"max_step unspecified and stop time unspecified. Setting max step to "
<<final_snapshot_fill_iteration<< " to fill all BTD snapshots." << std::endl;
warpx.updateMaxStep(final_snapshot_fill_iteration);
}
} else if (final_snapshot_fill_iteration > warpx.maxStep() || final_snapshot_fill_time > warpx.stopTime()) {
std::stringstream warn_string;
warn_string << "\nSimulation might not run long enough to fill all BTD snapshots.\n"
<< "Final step: " << warpx.maxStep() << "\n"
<<"Stop time: " << warpx.stopTime() << "\n"
<<"Last BTD snapshot fills around step: " << final_snapshot_fill_iteration << "\n"
<<" or time: " << final_snapshot_fill_time << "\n";
ablastr::warn_manager::WMRecordWarning(
"BTD", warn_string.str(),
ablastr::warn_manager::WarnPriority::low);
}
#ifdef WARPX_DIM_RZ
UpdateVarnamesForRZopenPMD();
#endif
}
void
BTDiagnostics::ReadParameters ()
{
BaseReadParameters();
WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::gamma_boost > 1.0_rt,
"gamma_boost must be > 1 to use the back-transformed diagnostics");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::boost_direction[2] == 1,
"The back transformed diagnostics currently only works if the boost is in the z-direction");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::do_moving_window,
"The moving window should be on if using the boosted frame diagnostic.");
// The next two asserts could be relaxed with respect to check to current step
WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::end_moving_window_step < 0,
"The moving window must not stop when using the boosted frame diagnostic.");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::start_moving_window_step == 0,
"The moving window must start at step zero for the boosted frame diagnostic.");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::moving_window_dir == WARPX_ZINDEX,
"The boosted frame diagnostic currently only works if the moving window is in the z direction.");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
m_format == "plotfile" || m_format == "openpmd",
"<diag>.format must be plotfile or openpmd for back transformed diagnostics");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
m_crse_ratio == amrex::IntVect(1),
"Only support for coarsening ratio of 1 in all directions is included for BTD\n"
);
// Read list of back-transform diag parameters requested by the user //
const amrex::ParmParse pp_diag_name(m_diag_name);
m_file_prefix = "diags/" + m_diag_name;
pp_diag_name.query("file_prefix", m_file_prefix);
pp_diag_name.query("do_back_transformed_fields", m_do_back_transformed_fields);
pp_diag_name.query("do_back_transformed_particles", m_do_back_transformed_particles);
AMREX_ALWAYS_ASSERT(m_do_back_transformed_fields or m_do_back_transformed_particles);
if (!m_do_back_transformed_fields) { m_varnames.clear(); }
std::vector<std::string> intervals_string_vec = {"0"};
bool const num_snapshots_specified = utils::parser::queryWithParser(
pp_diag_name, "num_snapshots_lab", m_num_snapshots_lab);
bool const intervals_specified = pp_diag_name.queryarr("intervals", intervals_string_vec);
if (num_snapshots_specified)
{
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!intervals_specified,
"For back-transformed diagnostics, user should specify either num_snapshots_lab or intervals, not both");
intervals_string_vec = {":" + std::to_string(m_num_snapshots_lab-1)};
}
m_intervals = utils::parser::BTDIntervalsParser(intervals_string_vec);
m_num_buffers = m_intervals.NumSnapshots();
// Read either dz_snapshots_lab or dt_snapshots_lab
bool snapshot_interval_is_specified = utils::parser::queryWithParser(
pp_diag_name, "dt_snapshots_lab", m_dt_snapshots_lab);
if ( utils::parser::queryWithParser(pp_diag_name, "dz_snapshots_lab", m_dz_snapshots_lab) ) {
m_dt_snapshots_lab = m_dz_snapshots_lab/PhysConst::c;
snapshot_interval_is_specified = true;
}
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(snapshot_interval_is_specified,
"For back-transformed diagnostics, user should specify either dz_snapshots_lab or dt_snapshots_lab");
utils::parser::queryWithParser(pp_diag_name, "buffer_size", m_buffer_size);
#ifdef WARPX_DIM_RZ
const amrex::Vector< std::string > BTD_varnames_supported = {"Er", "Et", "Ez",
"Br", "Bt", "Bz",
"jr", "jt", "jz", "rho"};
#else
const amrex::Vector< std::string > BTD_varnames_supported = {"Ex", "Ey", "Ez",
"Bx", "By", "Bz",
"jx", "jy", "jz", "rho"};
#endif
for (const auto& var : m_varnames) {
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
(utils::algorithms::is_in(BTD_varnames_supported, var )),
"Input error: field variable " + var + " in " + m_diag_name
+ ".fields_to_plot is not supported for BackTransformed diagnostics."
+ " Currently supported field variables for BackTransformed diagnostics "
+ "include Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, and rho in Cartesian coordinates and "
+ "Er, Et, Ez, Br, Bt, Bz, jr, jt, jz, and rho in cylindrical (RZ coordinates)");
}
const bool particle_fields_to_plot_specified = pp_diag_name.queryarr("particle_fields_to_plot", m_pfield_varnames);
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!particle_fields_to_plot_specified, "particle_fields_to_plot is currently not supported for BackTransformed Diagnostics");
if (m_varnames.empty()) {
m_do_back_transformed_fields = false;
}
}
bool
BTDiagnostics::DoDump (int step, int i_buffer, bool force_flush)
{
// Do not call dump if timestep < 0, i.e., at initialization time when step == -1
// or if the snapshot is already full and the files are closed.
if (step >= 0 && (m_snapshot_full[i_buffer] != 1)){
// If buffer for this lab snapshot is full then dump it and continue to collect
// slices afterwards
const auto is_buffer_full = buffer_full(i_buffer);
// or
// If last z-slice in the lab-frame snapshot is filled, call dump to
// write the buffer and close the file.
const auto last_z_slice_filled = (m_lastValidZSlice[i_buffer] == 1);
// or
// Do a forced dump at the end of the simulation, unless lab snapshot
// was already fully written and buffer was reset to zero size or that
// lab snapshot was not even started to be backtransformed yet
const auto do_forced_flush = (force_flush && !buffer_empty(i_buffer));
return is_buffer_full || last_z_slice_filled || do_forced_flush;
}
return false;
}
bool
BTDiagnostics::DoComputeAndPack (int step, bool force_flush)
{
// always set to true for BTDiagnostics since back-transform buffers are potentially
// computed and packed every timstep, except at initialization when step == -1, or when
// force_flush is set to true, because we dont need to redundantly re-compute
// buffers when force_flush = true. We only need to dump the buffers when
// force_flush=true. Note that the BTD computation is performed every timestep (step>=0)
return ( (step >=0 ) && (!force_flush) );
}
void
BTDiagnostics::InitializeBufferData ( int i_buffer , int lev, bool restart)
{
auto & warpx = WarpX::GetInstance();
// When restarting boosted simulations, the code below needs to take
// into account the fact that the position of the box at the beginning
// of the simulation, is not the one that we had at t=0 (because of the moving window)
const amrex::Real boosted_moving_window_v = (WarpX::moving_window_v - m_beta_boost*PhysConst::c)
/ (1._rt - m_beta_boost * WarpX::moving_window_v/PhysConst::c);
// Lab-frame time for the i^th snapshot
if (!restart) {
const amrex::Real zmax_0 = warpx.Geom(lev).ProbHi(m_moving_window_dir);
m_t_lab.at(i_buffer) = m_intervals.GetBTDIteration(i_buffer) * m_dt_snapshots_lab
+ m_gamma_boost*m_beta_boost*zmax_0/PhysConst::c;
}
// Define buffer domain in boosted frame at level, lev, with user-defined lo and hi
amrex::RealBox diag_dom;
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim ) {
// Setting lo-coordinate for the diag domain by taking the max of user-defined
// lo-cordinate and lo-coordinate of the simulation domain at level, lev
diag_dom.setLo(idim, std::max(m_lo[idim],warpx.Geom(lev).ProbLo(idim)) );
// Setting hi-coordinate for the diag domain by taking the max of user-defined
// hi-cordinate and hi-coordinate of the simulation domain at level, lev
diag_dom.setHi(idim, std::min(m_hi[idim],warpx.Geom(lev).ProbHi(idim)) );
}
// Initializing the m_buffer_box for the i^th snapshot.
// At initialization, the Box has the same index space as the boosted-frame
// As time-progresses, the z-dimension indices will be modified based on
// current_z_lab
amrex::IntVect lo(0);
amrex::IntVect hi(1);
for (int idim=0; idim < AMREX_SPACEDIM; ++idim) {
// lo index with same cell-size as simulation at level, lev.
const int lo_index = static_cast<int>( std::floor(
( diag_dom.lo(idim) - warpx.Geom(lev).ProbLo(idim) ) /
warpx.Geom(lev).CellSize(idim) ) );
// Taking max of (0,lo_index) because lo_index must always be >=0
lo[idim] = std::max( 0, lo_index );
// hi index with same cell-size as simulation at level, lev.
const int hi_index = static_cast<int>( std::ceil(
( diag_dom.hi(idim) - warpx.Geom(lev).ProbLo(idim) ) /
warpx.Geom(lev).CellSize(idim) ) );
// Taking max of (0,hi_index) because hi_index must always be >=0
// Subtracting by 1 because lo,hi indices are set to cell-centered staggering.
hi[idim] = std::max( 0, hi_index) - 1;
// if hi<=lo, then hi = lo + 1, to ensure one cell in that dimension
if ( hi[idim] <= lo[idim] ) {
hi[idim] = lo[idim] + 1;
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
m_crse_ratio[idim]==1, "coarsening ratio in reduced dimension must be 1."
);
}
}
const amrex::Box diag_box( lo, hi );
m_buffer_box[i_buffer] = diag_box;
m_snapshot_box[i_buffer] = diag_box;
// Define box array
amrex::BoxArray diag_ba(diag_box);
diag_ba.maxSize( warpx.maxGridSize( lev ) );
// Update the physical co-ordinates m_lo and m_hi using the final index values
// from the coarsenable, cell-centered BoxArray, ba.
for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
diag_dom.setLo( idim, warpx.Geom(lev).ProbLo(idim) +
diag_ba.getCellCenteredBox(0).smallEnd(idim) * warpx.Geom(lev).CellSize(idim));
diag_dom.setHi( idim, warpx.Geom(lev).ProbLo(idim) +
(diag_ba.getCellCenteredBox( static_cast<int>(diag_ba.size()-1) ).bigEnd(idim) + 1) * warpx.Geom(lev).CellSize(idim));
}
// Define buffer_domain in lab-frame for the i^th snapshot.
// Replace z-dimension with lab-frame co-ordinates.
const amrex::Real zmin_buffer_lab = ( diag_dom.lo(m_moving_window_dir) - boosted_moving_window_v * warpx.gett_new(0) )
/ ( (1.0_rt + m_beta_boost) * m_gamma_boost);
const amrex::Real zmax_buffer_lab = ( diag_dom.hi(m_moving_window_dir) - boosted_moving_window_v * warpx.gett_new(0) )
/ ( (1.0_rt + m_beta_boost) * m_gamma_boost);
// Initialize buffer counter and z-positions of the i^th snapshot in
// boosted-frame and lab-frame
m_buffer_counter[i_buffer] = 0;
m_current_z_lab[i_buffer] = 0._rt;
m_current_z_boost[i_buffer] = 0._rt;
// store old z boost before updated zboost position
m_old_z_boost[i_buffer] = m_current_z_boost[i_buffer];
// Now Update Current Z Positions
m_current_z_boost[i_buffer] = UpdateCurrentZBoostCoordinate(m_t_lab[i_buffer],
warpx.gett_new(lev) );
m_current_z_lab[i_buffer] = UpdateCurrentZLabCoordinate(m_t_lab[i_buffer],
warpx.gett_new(lev) );
// Compute number of cells in lab-frame required for writing Header file
// and potentially to generate Back-Transform geometry to ensure
// compatibility with plotfiles.
// For the z-dimension, number of cells in the lab-frame is
// computed using the coarsened cell-size in the lab-frame obtained using
// the ref_ratio at level, lev-1.
auto ref_ratio = amrex::IntVect(1);
if (lev > 0 ) { ref_ratio = WarpX::RefRatio(lev-1); }
// Number of lab-frame cells in z-direction at level, lev
const int num_zcells_lab = static_cast<int>( std::floor (
( zmax_buffer_lab - zmin_buffer_lab)
/ dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir])));
// Take the max of 0 and num_zcells_lab
const int Nz_lab = std::max( 0, num_zcells_lab );
#if (AMREX_SPACEDIM >= 2)
// Number of lab-frame cells in x-direction at level, lev
const int num_xcells_lab = static_cast<int>( std::floor (
( diag_dom.hi(0) - diag_dom.lo(0) )
/ warpx.Geom(lev).CellSize(0)
) );
// Take the max of 0 and num_ycells_lab
const int Nx_lab = std::max( 0, num_xcells_lab);
#endif
#if defined(WARPX_DIM_3D)
// Number of lab-frame cells in the y-direction at level, lev
const int num_ycells_lab = static_cast<int>( std::floor (
( diag_dom.hi(1) - diag_dom.lo(1) )
/ warpx.Geom(lev).CellSize(1)
) );
// Take the max of 0 and num_xcells_lab
const int Ny_lab = std::max( 0, num_ycells_lab );
m_snapshot_ncells_lab[i_buffer] = {Nx_lab, Ny_lab, Nz_lab};
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
m_snapshot_ncells_lab[i_buffer] = {Nx_lab, Nz_lab};
#else
m_snapshot_ncells_lab[i_buffer] = amrex::IntVect(Nz_lab);
#endif
// Box covering the extent of the user-defined diag in the back-transformed frame
// for the ith snapshot
// estimating the maximum number of buffer multifabs needed to obtain the
// full lab-frame snapshot
m_max_buffer_multifabs[i_buffer] = static_cast<int>( std::ceil (
amrex::Real(m_snapshot_ncells_lab[i_buffer][m_moving_window_dir]) /
amrex::Real(m_buffer_size) ) );
// number of cells in z is modified since each buffer multifab always
// contains a minimum m_buffer_size=256 cells
const int num_z_cells_in_snapshot = m_max_buffer_multifabs[i_buffer] * m_buffer_size;
if (!restart) {
m_snapshot_domain_lab[i_buffer] = diag_dom;
m_snapshot_domain_lab[i_buffer].setLo(m_moving_window_dir,
zmin_buffer_lab + WarpX::moving_window_v * m_t_lab[i_buffer]);
m_snapshot_domain_lab[i_buffer].setHi(m_moving_window_dir,
zmax_buffer_lab + WarpX::moving_window_v * m_t_lab[i_buffer]);
// To prevent round off errors, moving the snapshot domain by half a cell so that all the slices
// lie close to the cell-centers in the lab-frame grid instead of on the edge of cell.
const amrex::Real new_hi = m_snapshot_domain_lab[i_buffer].hi(m_moving_window_dir)
+ 0.5_rt * dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir]);
m_snapshot_domain_lab[i_buffer].setHi(m_moving_window_dir,new_hi);
const amrex::Real new_lo = m_snapshot_domain_lab[i_buffer].hi(m_moving_window_dir) -
num_z_cells_in_snapshot *
dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir]);
m_snapshot_domain_lab[i_buffer].setLo(m_moving_window_dir, new_lo);
}
// cell-centered index that corresponds to the hi-end of the lab-frame in the z-direction
// Adding 0.5 dz_lab so that we obtain the cell-centered index consistent to the hi-end
const int snapshot_kindex_hi = static_cast<int>(floor(
( m_snapshot_domain_lab[i_buffer].hi(m_moving_window_dir)
- (m_snapshot_domain_lab[i_buffer].lo(m_moving_window_dir)
+ 0.5*dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir])
)
) / dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir]) ));
m_snapshot_box[i_buffer].setBig( m_moving_window_dir, snapshot_kindex_hi);
m_snapshot_box[i_buffer].setSmall( m_moving_window_dir,
snapshot_kindex_hi - (num_z_cells_in_snapshot-1) );
// Setting hi k-index for the first buffer
if (!restart) {
m_buffer_k_index_hi[i_buffer] = m_snapshot_box[i_buffer].bigEnd(m_moving_window_dir);
}
}
void
BTDiagnostics::DefineCellCenteredMultiFab(int lev)
{
if (!m_do_back_transformed_fields) { return; }
// Creating MultiFab to store cell-centered data in boosted-frame for the entire-domain
// This MultiFab will store all the user-requested fields in the boosted-frame
auto & warpx = WarpX::GetInstance();
// The BoxArray is coarsened based on the user-defined coarsening ratio
amrex::BoxArray ba = warpx.boxArray(lev);
ba.coarsen(m_crse_ratio);
const amrex::DistributionMapping dmap = warpx.DistributionMap(lev);
const int ngrow = 1;
#ifdef WARPX_DIM_RZ
const int ncomps = WarpX::ncomps * static_cast<int>(m_cellcenter_varnames.size());
#else
const int ncomps = static_cast<int>(m_cellcenter_varnames.size());
#endif
WarpX::AllocInitMultiFab(m_cell_centered_data[lev], ba, dmap, ncomps, amrex::IntVect(ngrow), lev, "cellcentered_BTD", 0._rt);
}
void
BTDiagnostics::InitializeFieldFunctors (int lev)
{
// Initialize fields functors only if do_back_transformed_fields is selected
if (!m_do_back_transformed_fields) { return; }
#ifdef WARPX_DIM_RZ
// For RZ, initialize field functors RZ for openpmd
// This is a specialized call for intializing cell-center functors
// such that, all modes of a field component are stored contiguously
// For example, Er0, Er1_real, Er1_imag, etc
InitializeFieldFunctorsRZopenPMD(lev);
#else
auto & warpx = WarpX::GetInstance();
// Clear any pre-existing vector to release stored data
// This ensures that when domain is load-balanced, the functors point
// to the correct field-data pointers
m_all_field_functors[lev].clear();
// For back-transformed data, all the components are cell-centered and stored
// in a single multifab, m_cell_centered_data.
// Therefore, size of functors at all levels is 1.
const int num_BT_functors = 1;
m_all_field_functors[lev].resize(num_BT_functors);
m_cell_center_functors[lev].clear();
m_cell_center_functors[lev].resize( m_cellcenter_varnames.size() );
// Create an object of class BackTransformFunctor
for (int i = 0; i < num_BT_functors; ++i)
{
// coarsening ratio is not provided since the source MultiFab, m_cell_centered_data
// is coarsened based on the user-defined m_crse_ratio
const int nvars = static_cast<int>(m_varnames.size());
m_all_field_functors[lev][i] = std::make_unique<BackTransformFunctor>(
m_cell_centered_data[lev].get(), lev,
nvars, m_num_buffers, m_varnames, m_varnames_fields);
}
// Define all cell-centered functors required to compute cell-centere data
// Fill vector of cell-center functors for all field-components, namely,
// Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, and rho are included in the
// cell-center functors for BackTransform Diags
const auto m_cell_center_functors_at_lev_size = static_cast<int>(
m_cell_center_functors.at(lev).size());
for (int comp=0; comp<m_cell_center_functors_at_lev_size; comp++){
if ( m_cellcenter_varnames[comp] == "Ex" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Efield_aux(lev, 0), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "Ey" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Efield_aux(lev, 1), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "Ez" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Efield_aux(lev, 2), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "Bx" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 0), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "By" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 1), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "Bz" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 2), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "jx" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 0), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "jy" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 1), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "jz" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 2), lev, m_crse_ratio);
} else if ( m_cellcenter_varnames[comp] == "rho" ){
m_cell_center_functors[lev][comp] = std::make_unique<RhoFunctor>(lev, m_crse_ratio);
}
}
#endif
}
void
BTDiagnostics::UpdateVarnamesForRZopenPMD ()
{
#ifdef WARPX_DIM_RZ
auto & warpx = WarpX::GetInstance();
const int ncomp_multimodefab = warpx.get_pointer_Efield_aux(0,0)->nComp();
const int ncomp = ncomp_multimodefab;
const bool update_varnames = true;
if (update_varnames) {
const auto n_rz = ncomp * static_cast<int>(m_varnames_fields.size());
m_varnames.clear();
m_varnames.reserve(n_rz);
}
// AddRZ modes to output names for the back-transformed data
if (update_varnames) {
const auto m_varnames_fields_size = static_cast<int>(m_varnames_fields.size());
for (int comp=0; comp<m_varnames_fields_size; comp++)
{
if (m_varnames_fields[comp] == "Er") { AddRZModesToOutputNames(std::string("Er"), ncomp, false); }
if (m_varnames_fields[comp] == "Et") { AddRZModesToOutputNames(std::string("Et"), ncomp, false); }
if (m_varnames_fields[comp] == "Ez") { AddRZModesToOutputNames(std::string("Ez"), ncomp, false); }
if (m_varnames_fields[comp] == "Br") { AddRZModesToOutputNames(std::string("Br"), ncomp, false); }
if (m_varnames_fields[comp] == "Bt") { AddRZModesToOutputNames(std::string("Bt"), ncomp, false); }
if (m_varnames_fields[comp] == "Bz") { AddRZModesToOutputNames(std::string("Bz"), ncomp, false); }
if (m_varnames_fields[comp] == "jr") { AddRZModesToOutputNames(std::string("jr"), ncomp, false); }
if (m_varnames_fields[comp] == "jt") { AddRZModesToOutputNames(std::string("jt"), ncomp, false); }
if (m_varnames_fields[comp] == "jz") { AddRZModesToOutputNames(std::string("jz"), ncomp, false); }
if (m_varnames_fields[comp] == "rho") { AddRZModesToOutputNames(std::string("rho"),ncomp, false); }
}
}
// This function may be called multiple times, for different values of `lev`
// but the `varnames` need only be updated once.
const bool update_cellcenter_varnames = true;
if (update_cellcenter_varnames) {
const auto n_rz = ncomp * static_cast<int>(m_cellcenter_varnames.size());
m_cellcenter_varnames.clear();
m_cellcenter_varnames.reserve(n_rz);
const auto m_cellcenter_varnames_fields_size = static_cast<int>(m_cellcenter_varnames_fields.size());
for (int comp=0; comp<m_cellcenter_varnames_fields_size; comp++)
{
if ( m_cellcenter_varnames_fields[comp] == "Er" ) { AddRZModesToOutputNames(std::string("Er"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "Et" ) { AddRZModesToOutputNames(std::string("Et"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "Ez" ) { AddRZModesToOutputNames(std::string("Ez"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "Br" ) { AddRZModesToOutputNames(std::string("Br"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "Bt" ) { AddRZModesToOutputNames(std::string("Bt"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "Bz" ) { AddRZModesToOutputNames(std::string("Bz"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "jr" ) { AddRZModesToOutputNames(std::string("jr"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "jt" ) { AddRZModesToOutputNames(std::string("jt"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "jz" ) { AddRZModesToOutputNames(std::string("jz"), ncomp, true); }
if ( m_cellcenter_varnames_fields[comp] == "rho" ) { AddRZModesToOutputNames(std::string("rho"), ncomp, true); }
}
}
#endif
}
void
BTDiagnostics::InitializeFieldFunctorsRZopenPMD (int lev)
{
#ifdef WARPX_DIM_RZ
auto & warpx = WarpX::GetInstance();
const int ncomp_multimodefab = warpx.get_pointer_Efield_aux(0,0)->nComp();
const int ncomp = ncomp_multimodefab;
// Clear any pre-existing vector to release stored data
// This ensures that when domain is load-balanced, the functors point
// to the correct field-data pointers
m_all_field_functors[lev].clear();
// For back-transformed data, all the components are cell-centered and stored
// in a single multifab, m_cell_centered_data.
// Therefore, size of functors at all levels is 1
const int num_BT_functors = 1;
m_all_field_functors[lev].resize(num_BT_functors);
for (int i = 0; i < num_BT_functors; ++i) {
const int nvars = static_cast<int>(m_varnames.size());
m_all_field_functors[lev][i] = std::make_unique<BackTransformFunctor>(
m_cell_centered_data[lev].get(), lev,
nvars, m_num_buffers, m_varnames,
m_varnames_fields);
}
// Reset field functors for cell-center multifab
m_cell_center_functors[lev].clear();
m_cell_center_functors[lev].resize(m_cellcenter_varnames_fields.size());
const auto m_cell_center_functors_at_lev_size = static_cast<int>(m_cell_center_functors.at(lev).size());
for (int comp=0; comp<m_cell_center_functors_at_lev_size; comp++){
if ( m_cellcenter_varnames_fields[comp] == "Er" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Efield_aux(lev, 0), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "Et" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Efield_aux(lev, 1), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "Ez" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Efield_aux(lev, 2), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "Br" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 0), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "Bt" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 1), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "Bz" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 2), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "jr" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 0), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "jt" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 1), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "jz" ){
m_cell_center_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 2), lev, m_crse_ratio, false, ncomp);
} else if ( m_cellcenter_varnames_fields[comp] == "rho" ){
m_cell_center_functors[lev][comp] = std::make_unique<RhoFunctor>(lev, m_crse_ratio, -1, false, ncomp);
}
}
#endif
amrex::ignore_unused(lev);
}
void
BTDiagnostics::AddRZModesToOutputNames (const std::string& field, const int ncomp, bool cellcenter_data)
{
#ifdef WARPX_DIM_RZ
// In cylindrical geometry, real and imag part of each mode are also
// dumped to file separately, so they need to be added to m_varnames
// we number modes from 0 to (nmodes-1);
// mode 0 is purely real, all higher modes are complex
int const nmodes = (ncomp+1)/2;
if (cellcenter_data) {
m_cellcenter_varnames.push_back( field + "_0_real" );
for (int ic=1; ic < nmodes; ic++) {
m_cellcenter_varnames.push_back( field + "_" + std::to_string(ic) + "_real" );
m_cellcenter_varnames.push_back( field + "_" + std::to_string(ic) + "_imag" );
}
} else {
m_varnames.push_back(field + "_0_real");
for (int ic=1; ic < nmodes; ic++) {
m_varnames.push_back( field + "_" + std::to_string(ic) + "_real" );
m_varnames.push_back( field + "_" + std::to_string(ic) + "_imag" );
}
}
#else
amrex::ignore_unused(field, ncomp, cellcenter_data);
#endif
}
void
BTDiagnostics::PrepareBufferData ()
{
auto & warpx = WarpX::GetInstance();
const int num_BT_functors = 1;
for (int lev = 0; lev < nlev_output; ++lev)
{
for (int i = 0; i < num_BT_functors; ++i)
{
for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer )
{
m_old_z_boost[i_buffer] = m_current_z_boost[i_buffer];
// Update z-boost and z-lab positions
m_current_z_boost[i_buffer] = UpdateCurrentZBoostCoordinate(m_t_lab[i_buffer],
warpx.gett_new(lev) );
m_current_z_lab[i_buffer] = UpdateCurrentZLabCoordinate(m_t_lab[i_buffer],
warpx.gett_new(lev) );
}
}
}
}
void
BTDiagnostics::UpdateBufferData ()
{
const int num_BT_functors = 1;
for (int lev = 0; lev < nlev_output; ++lev)
{
for (int i = 0; i < num_BT_functors; ++i)
{
for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer )
{
const bool ZSliceInDomain = GetZSliceInDomainFlag (i_buffer, lev);
if (ZSliceInDomain) { ++m_buffer_counter[i_buffer]; }
// when the z-index is equal to the smallEnd of the snapshot box, then set lastValidZSlice to 1
if (k_index_zlab(i_buffer, lev) == m_snapshot_box[i_buffer].smallEnd(m_moving_window_dir)) {
m_lastValidZSlice[i_buffer] = 1;
}
}
}
}
}
void
BTDiagnostics::PrepareFieldDataForOutput ()
{
// Initialize fields functors only if do_back_transformed_fields is selected
if (!m_do_back_transformed_fields) { return; }
auto & warpx = WarpX::GetInstance();
// In this function, we will get cell-centered data for every level, lev,
// using the cell-center functors and their respective operators()
// Call m_cell_center_functors->operator
for (int lev = 0; lev < nmax_lev; ++lev) {
int icomp_dst = 0;
const auto m_cell_center_fuctors_at_lev_size = static_cast<int>(m_cell_center_functors.at(lev).size());
for (int icomp = 0; icomp<m_cell_center_fuctors_at_lev_size; ++icomp) {
// Call all the cell-center functors in m_cell_center_functors.
// Each of them computes cell-centered data for a field and
// stores it in cell-centered MultiFab, m_cell_centered_data[lev].
m_cell_center_functors[lev][icomp]->operator()(*m_cell_centered_data[lev], icomp_dst);
icomp_dst += m_cell_center_functors[lev][icomp]->nComp();
}
// Check that the proper number of user-requested components are cell-centered
AMREX_ALWAYS_ASSERT( icomp_dst == m_cellcenter_varnames.size() );
// fill boundary call is required to average_down (flatten) data to
// the coarsest level.
ablastr::utils::communication::FillBoundary(*m_cell_centered_data[lev], WarpX::do_single_precision_comms,
warpx.Geom(lev).periodicity());
}
// Flattening out MF over levels
for (int lev = warpx.finestLevel(); lev > 0; --lev) {
ablastr::coarsen::sample::Coarsen(*m_cell_centered_data[lev - 1], *m_cell_centered_data[lev], 0, 0,
static_cast<int>(m_cellcenter_varnames.size()), 0, WarpX::RefRatio(lev-1) );
}
const int num_BT_functors = 1;
for (int lev = 0; lev < nlev_output; ++lev)
{
for (int i = 0; i < num_BT_functors; ++i)
{
for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer )
{
// Check if the zslice is in domain
const bool ZSliceInDomain = GetZSliceInDomainFlag (i_buffer, lev);
// Initialize and define field buffer multifab if buffer is empty
const bool kindexInSnapshotBox = GetKIndexInSnapshotBoxFlag (i_buffer, lev);
if (kindexInSnapshotBox) {
if ( buffer_empty(i_buffer) ) {
if ( m_buffer_flush_counter[i_buffer] == 0 || m_first_flush_after_restart[i_buffer] == 1) {
// Compute the geometry, snapshot lab-domain extent
// and box-indices
DefineSnapshotGeometry(i_buffer, lev);
}
DefineFieldBufferMultiFab(i_buffer, lev);
}
}
if (ZSliceInDomain) {
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
m_current_z_lab[i_buffer] >= m_buffer_domain_lab[i_buffer].lo(m_moving_window_dir) and
m_current_z_lab[i_buffer] <= m_buffer_domain_lab[i_buffer].hi(m_moving_window_dir),
"z-slice in lab-frame (" +
std::to_string(m_current_z_lab[i_buffer]) +
") is outside the buffer domain physical extent (" +
std::to_string(m_buffer_domain_lab[i_buffer].lo(m_moving_window_dir)) +
" to " +
std::to_string(m_buffer_domain_lab[i_buffer].hi(m_moving_window_dir)) +
")."
);
}
m_all_field_functors[lev][i]->PrepareFunctorData (
i_buffer, ZSliceInDomain,
m_current_z_boost[i_buffer],
m_buffer_box[i_buffer],
k_index_zlab(i_buffer, lev),
m_snapshot_full[i_buffer] );
}
}
}
}
amrex::Real
BTDiagnostics::dz_lab (amrex::Real dt, amrex::Real ref_ratio) const
{
return PhysConst::c * dt * 1._rt/m_beta_boost * 1._rt/m_gamma_boost * 1._rt/ref_ratio;
}
int
BTDiagnostics::k_index_zlab (int i_buffer, int lev) const
{
auto & warpx = WarpX::GetInstance();
const amrex::Real prob_domain_zmin_lab = m_snapshot_domain_lab[i_buffer].lo( m_moving_window_dir );
auto ref_ratio = amrex::IntVect(1);
if (lev > 0 ) { ref_ratio = WarpX::RefRatio(lev-1); }
const int k_lab = static_cast<int>(std::floor (
( m_current_z_lab[i_buffer]
- (prob_domain_zmin_lab ) )
/ dz_lab( warpx.getdt(lev), ref_ratio[m_moving_window_dir] )
) ) + m_snapshot_box[i_buffer].smallEnd(m_moving_window_dir);
return k_lab;
}
void
BTDiagnostics::SetSnapshotFullStatus (const int i_buffer)
{
if (m_snapshot_full[i_buffer] == 1) { return; }
// if the last valid z-index of the snapshot, which is 0, is filled, then
// set the snapshot full integer to 1
if (m_lastValidZSlice[i_buffer] == 1) { m_snapshot_full[i_buffer] = 1; }
}
void
BTDiagnostics::DefineFieldBufferMultiFab (const int i_buffer, const int lev)
{
if (m_field_buffer_multifab_defined[i_buffer] == 1) { return; }
auto & warpx = WarpX::GetInstance();
const int hi_k_lab = m_buffer_k_index_hi[i_buffer];
m_buffer_box[i_buffer].setSmall( m_moving_window_dir, hi_k_lab - m_buffer_size + 1);
m_buffer_box[i_buffer].setBig( m_moving_window_dir, hi_k_lab );
amrex::BoxArray buffer_ba( m_buffer_box[i_buffer] );
// Generate a new distribution map for the back-transformed buffer multifab
const amrex::DistributionMapping buffer_dmap(buffer_ba);
// Number of guard cells for the output buffer is zero.
// Unlike FullDiagnostics, "m_format == sensei" option is not included here.
const int ngrow = 0;
m_mf_output[i_buffer][lev] = amrex::MultiFab( buffer_ba, buffer_dmap,
static_cast<int>(m_varnames.size()), ngrow );
m_mf_output[i_buffer][lev].setVal(0.);
auto ref_ratio = amrex::IntVect(1);
if (lev > 0 ) { ref_ratio = WarpX::RefRatio(lev-1); }
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
const amrex::Real cellsize = (idim < WARPX_ZINDEX)?
warpx.Geom(lev).CellSize(idim):
dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir]);
const amrex::Real buffer_lo = m_snapshot_domain_lab[i_buffer].lo(idim)
+ ( buffer_ba.getCellCenteredBox(0).smallEnd(idim)
- m_snapshot_box[i_buffer].smallEnd(idim)
) * cellsize;
const amrex::Real buffer_hi = m_snapshot_domain_lab[i_buffer].lo(idim)
+ ( buffer_ba.getCellCenteredBox( static_cast<int>(buffer_ba.size()-1) ).bigEnd(idim)
- m_snapshot_box[i_buffer].smallEnd(idim)
+ 1 ) * cellsize;
m_buffer_domain_lab[i_buffer].setLo(idim, buffer_lo);
m_buffer_domain_lab[i_buffer].setHi(idim, buffer_hi);
}
// Define the geometry object at level, lev, for the ith buffer.
if (lev == 0) {
// The extent of the physical domain covered by the ith buffer mf, m_mf_output
// Default non-periodic geometry for diags
amrex::Vector<int> BTdiag_periodicity(AMREX_SPACEDIM, 0);
// Box covering the extent of the user-defined diag in the back-transformed frame
const amrex::Box domain = buffer_ba.minimalBox();
// define the geometry object for the ith buffer using Physical co-ordinates
// of m_buffer_domain_lab[i_buffer].
m_geom_output[i_buffer][lev].define( domain, &m_buffer_domain_lab[i_buffer],
amrex::CoordSys::cartesian,
BTdiag_periodicity.data() );
} else if (lev > 0 ) {
// Refine the geometry object defined at the previous level, lev-1
m_geom_output[i_buffer][lev] = amrex::refine( m_geom_output[i_buffer][lev-1],
WarpX::RefRatio(lev-1) );
}
m_field_buffer_multifab_defined[i_buffer] = 1;
WARPX_ALWAYS_ASSERT_WITH_MESSAGE( m_mf_output[i_buffer][lev].boxArray().size() == 1,
"BoxArray size must be 1 for back-transformed diagnostics multifab that stores buffers");
}
void
BTDiagnostics::DefineSnapshotGeometry (const int i_buffer, const int lev)
{
if (m_snapshot_geometry_defined[i_buffer] == 1) { return; }
if (lev == 0) {
// Default non-periodic geometry for diags
amrex::Vector<int> BTdiag_periodicity(AMREX_SPACEDIM, 0);
// Define the geometry object for the ith snapshot using Physical co-ordinates
// of m_snapshot_domain_lab[i_buffer], that corresponds to the full snapshot
// in the back-transformed frame
m_geom_snapshot[i_buffer][lev].define( m_snapshot_box[i_buffer],
&m_snapshot_domain_lab[i_buffer],
amrex::CoordSys::cartesian,
BTdiag_periodicity.data() );
} else if (lev > 0) {
// Refine the geometry object defined at the previous level, lev-1
m_geom_snapshot[i_buffer][lev] = amrex::refine( m_geom_snapshot[i_buffer][lev-1],
WarpX::RefRatio(lev-1) );
}
m_snapshot_geometry_defined[i_buffer] = 1;
}
bool
BTDiagnostics::GetZSliceInDomainFlag (const int i_buffer, const int lev)
{
auto & warpx = WarpX::GetInstance();
const amrex::RealBox& boost_domain = warpx.Geom(lev).ProbDomain();
const amrex::Real buffer_zmin_lab = m_snapshot_domain_lab[i_buffer].lo( m_moving_window_dir );
const amrex::Real buffer_zmax_lab = m_snapshot_domain_lab[i_buffer].hi( m_moving_window_dir );
const bool slice_not_in_domain =
( m_current_z_boost[i_buffer] <= boost_domain.lo(m_moving_window_dir) ) ||
( m_current_z_boost[i_buffer] >= boost_domain.hi(m_moving_window_dir) ) ||
( m_current_z_lab[i_buffer] <= buffer_zmin_lab ) ||
( m_current_z_lab[i_buffer] >= buffer_zmax_lab );
return !slice_not_in_domain;
}
bool
BTDiagnostics::GetKIndexInSnapshotBoxFlag (const int i_buffer, const int lev)
{
return (k_index_zlab(i_buffer, lev) >= m_snapshot_box[i_buffer].smallEnd(m_moving_window_dir) &&
k_index_zlab(i_buffer, lev) <= m_snapshot_box[i_buffer].bigEnd(m_moving_window_dir));