forked from fmihpc/vlasiator
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathspatial_cell.hpp
1861 lines (1624 loc) · 85.7 KB
/
spatial_cell.hpp
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
/*
* This file is part of Vlasiator.
* Copyright 2010-2016 Finnish Meteorological Institute
*
* For details of usage, see the COPYING file and read the "Rules of the Road"
* at http://vlasiator.fmi.fi/
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along
* with this program; if not, write to the Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
/*!
Spatial cell class for Vlasiator that supports a variable number of velocity blocks.
*/
#ifndef VLASIATOR_SPATIAL_CELL_HPP
#define VLASIATOR_SPATIAL_CELL_HPP
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iostream>
#include <mpi.h>
#include <limits>
#include <stdint.h>
#include <vector>
#include <array>
#include <unordered_map>
#include <set>
#include <map>
#include <phiprof.hpp>
#include <tuple>
#include "memoryallocation.h"
#include "common.h"
#include "parameters.h"
#include "definitions.h"
#ifndef AMR
#include "velocity_mesh_old.h"
#else
#include "velocity_mesh_amr.h"
#endif
#include "amr_refinement_criteria.h"
#include "velocity_blocks.h"
#include "velocity_block_container.h"
#include "logger.h"
extern Logger logFile;
#ifndef NDEBUG
#define DEBUG_SPATIAL_CELL
#endif
typedef Parameters P;
// size of velocity blocks in velocity cells
#define block_vx_length WID
#define block_vy_length WID
#define block_vz_length WID
//this is also defined in common.h as SIZE_VELBLOCK, we should remove either one
#define VELOCITY_BLOCK_LENGTH WID3
//#define N_NEIGHBOR_VELOCITY_BLOCKS 28
/*!
Used as an error from functions returning velocity cells or
as a cell that would be outside of the velocity block
*/
#define error_velocity_cell 0xFFFFFFFFu
/*!
Used as an error from functions returning velocity cell indices or
as an index that would be outside of the velocity block
*/
#define error_velocity_cell_index 0xFFFFFFFFu
namespace spatial_cell {
namespace Transfer {
const uint64_t NONE = 0;
const uint64_t CELL_PARAMETERS = (1<<0);
const uint64_t CELL_DERIVATIVES = (1<<1);
const uint64_t VEL_BLOCK_LIST_STAGE1 = (1<<2);
const uint64_t VEL_BLOCK_LIST_STAGE2 = (1<<3);
const uint64_t VEL_BLOCK_DATA = (1<<4);
const uint64_t VEL_BLOCK_PARAMETERS = (1<<6);
const uint64_t VEL_BLOCK_WITH_CONTENT_STAGE1 = (1<<7);
const uint64_t VEL_BLOCK_WITH_CONTENT_STAGE2 = (1<<8);
const uint64_t CELL_SYSBOUNDARYFLAG = (1<<9);
const uint64_t CELL_E = (1<<10);
const uint64_t CELL_EDT2 = (1<<11);
const uint64_t CELL_PERB = (1<<12);
const uint64_t CELL_PERBDT2 = (1<<13);
const uint64_t CELL_BGB = (1<<14);
const uint64_t CELL_RHO_RHOV = (1<<15);
const uint64_t CELL_RHODT2_RHOVDT2 = (1<<16);
const uint64_t CELL_BVOL = (1<<17);
const uint64_t CELL_BVOL_DERIVATIVES = (1<<18);
const uint64_t CELL_DIMENSIONS = (1<<19);
const uint64_t CELL_IOLOCALCELLID = (1<<20);
const uint64_t NEIGHBOR_VEL_BLOCK_DATA = (1<<21);
const uint64_t CELL_HALL_TERM = (1<<22);
const uint64_t CELL_P = (1<<23);
const uint64_t CELL_PDT2 = (1<<24);
const uint64_t CELL_RHOQ_TOT = (1<<25);
const uint64_t CELL_PHI = (1<<26);
const uint64_t POP_METADATA = (1<<27);
const uint64_t RANDOMGEN = (1<<28);
const uint64_t CELL_GRADPE_TERM = (1<<29);
//all data
const uint64_t ALL_DATA =
CELL_PARAMETERS
| CELL_DERIVATIVES | CELL_BVOL_DERIVATIVES
| VEL_BLOCK_DATA
| CELL_SYSBOUNDARYFLAG
| POP_METADATA | RANDOMGEN;
//all data, except the distribution function
const uint64_t ALL_SPATIAL_DATA =
CELL_PARAMETERS
| CELL_DERIVATIVES | CELL_BVOL_DERIVATIVES
| CELL_SYSBOUNDARYFLAG
| POP_METADATA | RANDOMGEN;
}
typedef std::array<unsigned int, 3> velocity_cell_indices_t; /**< Defines the indices of a velocity cell in a velocity block.
* Indices start from 0 and the first value is the index in x direction.
* Note: these are the (i,j,k) indices of the cell within the block.
* Valid values are ([0,block_vx_length[,[0,block_vy_length[,[0,block_vz_length[).*/
typedef std::array<vmesh::LocalID,3> velocity_block_indices_t; /**< Defines the indices of a velocity block in the velocity grid.
* Indices start from 0 and the first value is the index in x direction.
* Note: these are the (i,j,k) indices of the block.
* Valid values are ([0,vx_length[,[0,vy_length[,[0,vz_length[).*/
/** Wrapper for variables needed for each particle species.*/
struct Population {
Real max_dt[2]; /**< Element[0] is max_r_dt, element[1] max_v_dt.*/
vmesh::LocalID N_blocks; /**< Number of velocity blocks, used when receiving velocity
* mesh from remote neighbors using MPI.*/
Real velocityBlockMinValue;
vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID> vmesh; /**< Velocity mesh. Contains all velocity blocks that exist
* in this spatial cell. Cells are identified by their unique
* global IDs.*/
vmesh::VelocityBlockContainer<vmesh::LocalID> blockContainer; /**< Velocity block data.*/
};
class SpatialCell {
public:
SpatialCell();
SpatialCell(const SpatialCell& other);
// Following functions return velocity grid metadata //
template<int PAD> void fetch_data(const vmesh::GlobalID& blockGID,const vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& vmesh,
const Realf* src,Realf* array);
template<int PAD> void fetch_acc_data(const vmesh::GlobalID& blockGID,const int& dim,
vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& vmesh,
const Realf* src,Realf* array,Real cellSizeFractions[2]);
vmesh::GlobalID find_velocity_block(uint8_t& refLevel,vmesh::GlobalID cellIndices[3],const int& popID);
Realf* get_data(const int& popID);
const Realf* get_data(const int& popID) const;
Realf* get_data(const vmesh::LocalID& blockLID,const int& popID);
const Realf* get_data(const vmesh::LocalID& blockLID,const int& popID) const;
Real* get_block_parameters(const int& popID);
const Real* get_block_parameters(const int& popID) const;
Real* get_block_parameters(const vmesh::LocalID& blockLID,const int& popID);
const Real* get_block_parameters(const vmesh::LocalID& blockLID,const int& popID) const;
Real* get_cell_parameters();
const Real* get_cell_parameters() const;
vmesh::LocalID get_number_of_velocity_blocks(const int& popID) const;
vmesh::LocalID get_number_of_all_velocity_blocks() const;
int get_number_of_populations() const;
uint8_t get_maximum_refinement_level(const int& popID);
const Real& get_max_r_dt(const int& popID) const;
const Real& get_max_v_dt(const int& popID) const;
const vmesh::LocalID* get_velocity_grid_length(const int& popID,const uint8_t& refLevel=0);
const Real* get_velocity_grid_block_size(const int& popID,const uint8_t& refLevel=0);
const Real* get_velocity_grid_cell_size(const int& popID,const uint8_t& refLevel=0);
void get_velocity_block_coordinates(const int& popID,const vmesh::GlobalID& globalID,Real* coords);
velocity_block_indices_t get_velocity_block_indices(const int& popID,const vmesh::GlobalID globalID); // OK
velocity_block_indices_t get_velocity_block_indices(const int& popID,const vmesh::GlobalID globalID,uint8_t& refLevel);
vmesh::GlobalID get_velocity_block(const int& popID,vmesh::GlobalID blockIndices[3],const uint8_t& refLevel) const;
vmesh::GlobalID get_velocity_block(const int& popID,const velocity_block_indices_t indices,const uint8_t& refLevel) const;
vmesh::GlobalID get_velocity_block(const int& popID,const Real* coords,const uint8_t& refLevel=0) const;
vmesh::GlobalID get_velocity_block(const int& popID,const Real vx,const Real vy,const Real vz,const uint8_t& refLevel=0) const;
vmesh::GlobalID get_velocity_block_child(const int& popID,const vmesh::GlobalID& blockGID,const uint8_t& refLevel,
const int& i_cell,const int& j_cell,const int& k_cell);
void get_velocity_block_children_local_ids(const vmesh::GlobalID& blockGID,
std::vector<vmesh::LocalID>& childrenLIDs,
const int& popID);
vmesh::GlobalID get_velocity_block_parent(const int& popID,const vmesh::GlobalID& blockGID);
vmesh::GlobalID get_velocity_block_global_id(const vmesh::LocalID& blockLID,const int& popID) const;
vmesh::LocalID get_velocity_block_local_id(const vmesh::GlobalID& blockGID,const int& popID) const;
void get_velocity_block_size(const int& popID,const vmesh::GlobalID block,Real size[3]);
Real get_velocity_block_vx_min(const int& popID,const vmesh::GlobalID block) const;
Real get_velocity_block_vx_max(const int& popID,const vmesh::GlobalID block) const;
Real get_velocity_block_vy_min(const int& popID,const vmesh::GlobalID block) const;
Real get_velocity_block_vy_max(const int& popID,const vmesh::GlobalID block) const;
Real get_velocity_block_vz_min(const int& popID,const vmesh::GlobalID block) const;
Real get_velocity_block_vz_max(const int& popID,const vmesh::GlobalID block) const;
velocity_cell_indices_t get_velocity_cell_indices(const unsigned int cell) const;
unsigned int get_velocity_cell(const velocity_cell_indices_t indices) const;
unsigned int get_velocity_cell(const int& popID,const vmesh::GlobalID velocity_block,const Real vx,const Real vy,const Real vz) const;
Real get_velocity_cell_vx_min(const int& popID,const vmesh::GlobalID velocity_block,const unsigned int velocity_cell) const;
Real get_velocity_cell_vx_max(const int& popID,const vmesh::GlobalID velocity_block,const unsigned int velocity_cell) const;
Real get_velocity_cell_vy_min(const int& popID,const vmesh::GlobalID velocity_block,const unsigned int velocity_cell) const;
Real get_velocity_cell_vy_max(const int& popID,const vmesh::GlobalID velocity_block,const unsigned int velocity_cell) const;
Real get_velocity_cell_vz_min(const int& popID,const vmesh::GlobalID velocity_block,const unsigned int velocity_cell) const;
Real get_velocity_cell_vz_max(const int& popID,const vmesh::GlobalID velocity_block,const unsigned int velocity_cell) const;
const Real* get_velocity_grid_min_limits(const int& popID);
const Real* get_velocity_grid_max_limits(const int& popID);
bool initialize_mesh();
static unsigned int invalid_block_index();
static vmesh::GlobalID invalid_global_id();
static vmesh::LocalID invalid_local_id();
size_t count(const vmesh::GlobalID& block,const int& popID) const;
void add_values(const vmesh::GlobalID& targetGID,
std::unordered_map<vmesh::GlobalID,Realf[(WID+2)*(WID+2)*(WID+2)]>& sourceData,
const int& popID);
void printMeshSizes();
static bool setCommunicatedSpecies(const int& popID);
// Following functions adjust velocity blocks stored on the cell //
bool add_velocity_block(const vmesh::GlobalID& block,const int& popID);
void add_velocity_blocks(const std::vector<vmesh::GlobalID>& blocks,const int& popID);
bool add_velocity_block_octant(const vmesh::GlobalID& blockGID,const int& popID);
void adjustSingleCellVelocityBlocks(const int& popID);
void adjust_velocity_blocks(const std::vector<SpatialCell*>& spatial_neighbors,
const int& popID,
bool doDeleteEmptyBlocks=true);
void update_velocity_block_content_lists(const int& popID);
bool checkMesh(const int& popID);
void clear(const int& popID);
void coarsen_block(const vmesh::GlobalID& parent,const std::vector<vmesh::GlobalID>& children,const int& popID);
void coarsen_blocks(amr_ref_criteria::Base* evaluator,const int& popID);
uint64_t get_cell_memory_capacity();
uint64_t get_cell_memory_size();
void merge_values(const int& popID);
void prepare_to_receive_blocks(const int& popID);
bool shrink_to_fit();
size_t size(const int& popID) const;
void remove_velocity_block(const vmesh::GlobalID& block,const int& popID);
void swap(vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& vmesh,
vmesh::VelocityBlockContainer<vmesh::LocalID>& blockContainer,const int& popID);
vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& get_velocity_mesh(const size_t& popID);
vmesh::VelocityBlockContainer<vmesh::LocalID>& get_velocity_blocks(const size_t& popID);
vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& get_velocity_mesh_temporary();
vmesh::VelocityBlockContainer<vmesh::LocalID>& get_velocity_blocks_temporary();
Realf get_value(const Real vx,const Real vy,const Real vz,const int& popID) const;
Realf get_value(const vmesh::GlobalID& blockGID, const unsigned int cell, const int& popID) const;
void increment_value(const Real vx,const Real vy,const Real vz,const Realf value,const int& popID);
void increment_value(const vmesh::GlobalID& block,const unsigned int cell,const Realf value,const int& popID);
void set_max_r_dt(const int& popID,const Real& value);
void set_max_v_dt(const int& popID,const Real& value);
void set_value(const Real vx, const Real vy, const Real vz, const Realf value,const int& popID);
void set_value(const vmesh::GlobalID& block,const unsigned int cell, const Realf value,const int& popID);
void refine_block(const vmesh::GlobalID& block,std::map<vmesh::GlobalID,vmesh::LocalID>& insertedBlocks,
const int& popID);
bool velocity_block_has_children(const vmesh::GlobalID& blockGID,const int& popID) const;
vmesh::GlobalID velocity_block_has_grandparent(const vmesh::GlobalID& blockGID,const int& popID) const;
// Following functions are related to MPI //
std::tuple<void*, int, MPI_Datatype> get_mpi_datatype(const CellID cellID,const int sender_rank,const int receiver_rank,
const bool receiving,const int neighborhood);
static uint64_t get_mpi_transfer_type(void);
static void set_mpi_transfer_type(const uint64_t type,bool atSysBoundaries=false);
void set_mpi_transfer_enabled(bool transferEnabled);
void updateSparseMinValue(const int& popID);
Real getVelocityBlockMinValue(const int& popID) const;
// Random number generator functions
//char* get_rng_state_buffer();
//random_data* get_rng_data_buffer();
// Member variables //
Real derivatives[fieldsolver::N_SPATIAL_CELL_DERIVATIVES]; /**< Derivatives of bulk variables in this spatial cell.*/
Real derivativesBVOL[bvolderivatives::N_BVOL_DERIVATIVES]; /**< Derivatives of BVOL needed by the acceleration.
* Separate array because it does not need to be communicated.*/
Real parameters[CellParams::N_SPATIAL_CELL_PARAMS]; /**< Bulk variables in this spatial cell.*/
Realf null_block_data[WID3];
uint64_t ioLocalCellId; /**< Local cell ID used for IO, not needed elsewhere
* and thus not being kept up-to-date.*/
//vmesh::LocalID mpi_number_of_blocks; /**< Number of blocks in mpi_velocity_block_list.*/
Realf* neighbor_block_data; /**< Pointers for translation operator. We can point to neighbor
* cell block data. We do not allocate memory for the pointer.*/
vmesh::LocalID neighbor_number_of_blocks;
uint sysBoundaryFlag; /**< What type of system boundary does the cell belong to.
* Enumerated in the sysboundarytype namespace's enum.*/
uint sysBoundaryLayer; /**< Layers counted from closest systemBoundary. If 0 then it has not
* been computed. First sysboundary layer is layer 1.*/
int sysBoundaryLayerNew;
std::vector<vmesh::GlobalID> velocity_block_with_content_list; /**< List of existing cells with content, only up-to-date after
* call to update_has_content().*/
vmesh::LocalID velocity_block_with_content_list_size; /**< Size of vector. Needed for MPI communication of size before actual list transfer.*/
std::vector<vmesh::GlobalID> velocity_block_with_no_content_list; /**< List of existing cells with no content, only up-to-date after
* call to update_has_content. This is also never transferred
* over MPI, so is invalid on remote cells.*/
static uint64_t mpi_transfer_type; /**< Which data is transferred by the mpi datatype given by spatial cells.*/
static bool mpiTransferAtSysBoundaries; /**< Do we only transfer data at boundaries (true), or in the whole system (false).*/
private:
SpatialCell& operator=(const SpatialCell&);
bool compute_block_has_content(const vmesh::GlobalID& block,const int& popID) const;
void merge_values_recursive(const int& popID,vmesh::GlobalID parentGID,vmesh::GlobalID blockGID,uint8_t refLevel,bool recursive,const Realf* data,
std::set<vmesh::GlobalID>& blockRemovalList);
static int activePopID;
bool initialized;
bool mpiTransferEnabled;
// Random number generator state variables, used for running reproducible
// simulations that do not depend on the number of threads of MPI processes used.
//char rngStateBuffer[256]; /**< Random number generator state buffer.*/
//random_data rngDataBuffer; /**< Random number generator data buffer.*/
// Temporary mesh used in acceleration and propagation.
vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID> vmeshTemp; /**< Temporary velocity mesh that is used in Vlasov solver.
* NOTE: Do not call the get-functions using this mesh as object
* before you have set the correct meshID using setMesh function.*/
vmesh::VelocityBlockContainer<vmesh::LocalID> blockContainerTemp;
std::vector<spatial_cell::Population> populations; /**< Particle population variables.*/
};
/****************************
* Velocity block functions *
****************************/
template<int PAD> inline
void SpatialCell::fetch_acc_data(const vmesh::GlobalID& blockGID,const int& dim,
vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& vmesh,
const Realf* src,Realf* array,Real cellSizeFractions[2]) {
const vmesh::LocalID blockLID = vmesh.getLocalID(blockGID);
#ifdef DEBUG_SPATIAL_CELL
if (blockGID == vmesh.invalidGlobalID() || blockLID == vmesh.invalidLocalID()) {
std::cerr << "ERROR: block has invalid global or local index " << __FILE__ << ':' << __LINE__ << std::endl;
exit(1);
}
#endif
const Realf* ptr = NULL;
uint8_t refLevel;
vmesh::LocalID i_block,j_block,k_block;
vmesh.getIndices(blockGID,refLevel,i_block,j_block,k_block);
// Copy values from x face neighbors:
std::vector<vmesh::LocalID> nbrIDs;
int32_t refLevelDiff;
Real crd;
switch (dim) {
case 0: // Transpose i->k, j->j, k->i
ptr = src + blockLID*WID3; // Copy values from this block
for (int k=0; k<WID; ++k) for (int j=0; j<WID; ++j) for (int i=0; i<WID; ++i) {
array[vblock::index(k,j,i+PAD)] = ptr[vblock::index(i,j,k)];
}
for (int i_nbr_off=-1; i_nbr_off<2; i_nbr_off+=2) { // Copy values from x face neighbors:
// Get local IDs of neighbor blocks
vmesh.getNeighborsExistingAtOffset(blockGID,i_nbr_off,+0,+0,nbrIDs,refLevelDiff);
// Position that is used to interpolate values from neighbor blocks
Real pos[3];
if (i_nbr_off < 0) crd = WID-0.5-(PAD-1);
else crd = 0.5;
// i-index to array where interpolated values are stored
uint32_t i_trgt = 0;
if (i_nbr_off > 0) i_trgt = WID+PAD;
if (nbrIDs.size() > 0) { // This block has at least one existing neighbor
if (refLevelDiff == -1) { // Neighbor is one level coarser, interpolate
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data; // (this check might not be necessary here)
else ptr = src + nbrIDs[0]*WID3;
for (uint32_t i=0; i<PAD; ++i) for (uint32_t k=0; k<WID; ++k) for (uint32_t j=0; j<WID; ++j) {
pos[0] = crd + i;
pos[1] = 2*(j_block%2) + j/2 + 0.5;
pos[2] = 2*(k_block%2) + k/2 + 0.5;
array[vblock::index(k,j,i_trgt+i)] = vblock::interp_xy<vblock::interpmethod::NGP>(pos,ptr);
}
cellSizeFractions[(i_nbr_off+1)/2] = 2.0;
} else if (refLevelDiff == 0) { // Neighbor at same level, copy data
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data; // (this check might not be necessary here)
else ptr = src + nbrIDs[0]*WID3;
uint32_t i_src = 0;
if (i_nbr_off < 0) i_src = WID-PAD;
for (uint32_t i=0; i<PAD; ++i) for (uint32_t k=0; k<WID; ++k) for (uint32_t j=0; j<WID; ++j) {
array[vblock::index(k,j,i_trgt+i)] = ptr[vblock::index(i_src+i,j,k)];
}
cellSizeFractions[(i_nbr_off+1)/2] = 1.0;
} else if (refLevelDiff == +1) { // nbr one level more refined, interpolate from four neighbors
for (uint32_t i=0; i<PAD; ++i) for (uint32_t k=0; k<WID; ++k) for (uint32_t j=0; j<WID; ++j) {
int index = (k/2)*2 + j/2;
if (nbrIDs[index] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[index]*WID3;
pos[0] = crd + i;
pos[1] = 2*(j%2) + 1;
pos[2] = 2*(k%2) + 1;
array[vblock::index(k,j,i_trgt+i)] = vblock::interp_xy<vblock::interpmethod::CIC>(pos,ptr);
}
cellSizeFractions[(i_nbr_off+1)/2] = 0.5;
}
} else { // Neighbor does not exist, return zero values
for (uint32_t i=0; i<PAD; ++i) for (uint32_t k=0; k<WID; ++k) for (uint32_t j=0; j<WID; ++j) {
array[vblock::index(k,j,i_trgt+i)] = 0.0;
}
cellSizeFractions[(i_nbr_off+1)/2] = 1.0;
}
}
break;
case 1: // Transpose i->i, j->k, k->j
ptr = src + blockLID*WID3; // Copy values from this block
for (int k=0; k<WID; ++k) for (int j=0; j<WID; ++j) for (int i=0; i<WID; ++i) {
array[vblock::index(i,k,j+PAD)] = ptr[vblock::index(i,j,k)];
}
for (int j_nbr_off=-1; j_nbr_off<2; j_nbr_off+=2) { // Copy values from y face neighbors:
// Get local IDs of neighbor blocks
vmesh.getNeighborsExistingAtOffset(blockGID,+0,j_nbr_off,+0,nbrIDs,refLevelDiff);
// Position that is used to interpolate values from neighbor blocks
Real pos[3];
if (j_nbr_off < 0) crd = WID-0.5-(PAD-1);
else crd = 0.5;
// j-index to array where interpolated values are stored
uint32_t j_trgt = 0;
if (j_nbr_off > 0) j_trgt = WID+PAD;
if (nbrIDs.size() > 0) { // This block has at least one existing neighbor
if (refLevelDiff == -1) { // Neighbor is one level coarser, interpolate
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data; // (this check might not be necessary here)
else ptr = src + nbrIDs[0]*WID3;
for (uint32_t j=0; j<PAD; ++j) for (uint32_t k=0; k<WID; ++k) for (uint32_t i=0; i<WID; ++i) {
pos[0] = 2*(i_block%2) + i/2 + 0.5;
pos[1] = crd + j;
pos[2] = 2*(k_block%2) + k/2 + 0.5;
array[vblock::index(i,k,j_trgt+j)] = vblock::interp_xy<vblock::interpmethod::NGP>(pos,ptr);
}
cellSizeFractions[(j_nbr_off+1)/2] = 2.0;
} else if (refLevelDiff == 0) { // Neighbor at same level, copy data
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data; // (this check might not be necessary here)
else ptr = src + nbrIDs[0]*WID3;
uint32_t j_src = 0;
if (j_nbr_off < 0) j_src = WID-PAD;
for (uint32_t j=0; j<PAD; ++j) for (uint32_t k=0; k<WID; ++k) for (uint32_t i=0; i<WID; ++i) {
array[vblock::index(i,k,j_trgt+j)] = ptr[vblock::index(i,j_src+j,k)];
}
cellSizeFractions[(j_nbr_off+1)/2] = 1.0;
} else if (refLevelDiff == +1) { // nbr one level more refined, interpolate from four neighbors
for (uint32_t j=0; j<PAD; ++j) for (uint32_t k=0; k<WID; ++k) for (uint32_t i=0; i<WID; ++i) {
// Iterate over the four neighbors. If the neighbor does not exist,
// interpolate values from the null block
int index = (k/2)*2 + i/2;
if (nbrIDs[index] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[index]*WID3;
pos[0] = 2*(i%2) + 1;
pos[1] = crd + j;
pos[2] = 2*(k%2) + 1;
array[vblock::index(i,k,j_trgt+j)] = vblock::interp_xy<vblock::interpmethod::CIC>(pos,ptr);
}
cellSizeFractions[(j_nbr_off+1)/2] = 0.5;
}
} else { // Neighbor does not exist, return zero values
for (uint32_t j=0; j<PAD; ++j) for (uint32_t k=0; k<WID; ++k) for (uint32_t i=0; i<WID; ++i) {
array[vblock::index(i,k,j_trgt+j)] = 0.0;
}
cellSizeFractions[(j_nbr_off+1)/2] = 1.0;
}
}
break;
case 2:
ptr = src + blockLID*WID3; // Copy values from this block
for (int k=0; k<WID; ++k) for (int j=0; j<WID; ++j) for (int i=0; i<WID; ++i) {
array[vblock::index(i,j,k+PAD)] = ptr[vblock::index(i,j,k)];
}
for (int k_nbr_off=-1; k_nbr_off<2; k_nbr_off+=2) { // Copy values from z face neighbors:
// Get local IDs of neighbor blocks
vmesh.getNeighborsExistingAtOffset(blockGID,+0,+0,k_nbr_off,nbrIDs,refLevelDiff);
// Position that is used to interpolate values from neighbor blocks
Real pos[3];
if (k_nbr_off < 0) crd = WID-0.5-(PAD-1);
else crd = 0.5;
// k-index to array where interpolated values are stored
uint32_t k_trgt = 0;
if (k_nbr_off > 0) k_trgt = WID+PAD;
if (nbrIDs.size() > 0) { // This block has at least one existing neighbor
if (refLevelDiff == -1) { // Neighbor is one level coarser, interpolate
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data; // (this check might not be necessary here)
else ptr = src + nbrIDs[0]*WID3;
for (uint32_t k=0; k<PAD; ++k) for (uint32_t j=0; j<WID; ++j) for (uint32_t i=0; i<WID; ++i) {
pos[0] = 2*(i_block%2) + i/2 + 0.5;
pos[1] = 2*(j_block%2) + j/2 + 0.5;
pos[2] = crd + k;
array[vblock::index(i,j,k_trgt+k)] = vblock::interp_xy<vblock::interpmethod::NGP>(pos,ptr);
}
cellSizeFractions[(k_nbr_off+1)/2] = 2.0;
} else if (refLevelDiff == 0) { // Neighbor at same level, copy data
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data; // (this check might not be necessary here)
else ptr = src + nbrIDs[0]*WID3;
uint32_t k_src = 0;
if (k_nbr_off < 0) k_src = WID-PAD;
for (uint32_t k=0; k<PAD; ++k) for (uint32_t j=0; j<WID; ++j) for (uint32_t i=0; i<WID; ++i) {
array[vblock::index(i,j,k_trgt+k)] = ptr[vblock::index(i,j,k_src+k)];
}
cellSizeFractions[(k_nbr_off+1)/2] = 1.0;
} else if (refLevelDiff == +1) { // nbr one level more refined, interpolate from four neighbors
for (uint32_t k=0; k<PAD; ++k) for (uint32_t j=0; j<WID; ++j) for (uint32_t i=0; i<WID; ++i) {
// Iterate over the four neighbors. If the neighbor does not exist,
// interpolate values from the null block
int index = (j/2)*2 + i/2;
if (nbrIDs[index] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[index]*WID3;
pos[0] = 2*(i%2) + 1;
pos[1] = 2*(j%2) + 1;
pos[2] = crd + k;
array[vblock::index(i,j,k_trgt+k)] = vblock::interp_xy<vblock::interpmethod::CIC>(pos,ptr);
}
cellSizeFractions[(k_nbr_off+1)/2] = 0.5;
}
} else { // Neighbor does not exist, return zero values
for (uint32_t k=0; k<PAD; ++k) for (uint32_t j=0; j<WID; ++j) for (uint32_t i=0; i<WID; ++i) {
array[vblock::index(i,j,k_trgt+k)] = 0.0;
}
cellSizeFractions[(k_nbr_off+1)/2] = 1.0;
}
}
break;
} // end switch
}
template<int PAD> inline
void SpatialCell::fetch_data(const vmesh::GlobalID& blockGID,
const vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>& vmesh,
const Realf* src,Realf* array) {
//const vmesh::LocalID blockLID = get_velocity_block_local_id(blockGID);
const vmesh::LocalID blockLID = vmesh.getLocalID(blockGID);
if (blockLID == invalid_local_id()) {
std::cerr << "ERROR: invalid local id in " << __FILE__ << ' ' << __LINE__ << std::endl;
exit(1);
}
// Copy values from this block:
const Realf* ptr = src + blockLID*WID3;
for (int k=0; k<WID; ++k) for (int j=0; j<WID; ++j) for (int i=0; i<WID; ++i) {
array[vblock::padIndex<PAD>(i+PAD,j+PAD,k+PAD)] = ptr[vblock::index(i,j,k)];
}
uint8_t refLevel;
vmesh::LocalID i_block,j_block,k_block;
vmesh.getIndices(blockGID,refLevel,i_block,j_block,k_block);
// Copy values from x face neighbors:
std::vector<vmesh::LocalID> nbrIDs;
int32_t refLevelDiff;
Real crd;
for (int i_nbr_off=-1; i_nbr_off<2; i_nbr_off+=2) {
vmesh.getNeighborsExistingAtOffset(blockGID,i_nbr_off,+0,+0,nbrIDs,refLevelDiff);
Real pos[3];
if (i_nbr_off < 0) crd = WID-0.5-(PAD-1);
else crd = 0.5;
uint32_t i_trgt = 0;
if (i_nbr_off > 0) i_trgt = WID+PAD;
if (nbrIDs.size() > 0) {
if (refLevelDiff == -1) { // nbr one level coarser, interpolate
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[0]*WID3;
for (uint32_t i=0; i<PAD; ++i) for (uint32_t k=0; k<WID; ++k) for (uint32_t j=0; j<WID; ++j) {
pos[0] = crd + i;
pos[1] = 2*(j_block%2) + j/2 + 0.5;
pos[2] = 2*(k_block%2) + k/2 + 0.5;
array[vblock::padIndex<PAD>(i_trgt+i,j+PAD,k+PAD)] = vblock::interp_yz<vblock::interpmethod::NGP>(pos,ptr);
}
} else if (refLevelDiff == 0) { // nbr at same level, simple data copy
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[0]*WID3;
uint32_t i_src = 0;
if (i_nbr_off < 0) i_src = WID-PAD;
for (uint32_t i=0; i<PAD; ++i) for (uint32_t k=0; k<WID; ++k) for (uint32_t j=0; j<WID; ++j) {
array[vblock::padIndex<PAD>(i_trgt+i,j+PAD,k+PAD)] = ptr[vblock::index(i_src+i,j,k)];
}
} else if (refLevelDiff == +1) { // nbr one level more refined, interpolate from four neighbors
for (uint32_t i=0; i<PAD; ++i) for (uint32_t k=0; k<WID; ++k) for (uint32_t j=0; j<WID; ++j) {
int index = (k/2)*2 + j/2;
if (nbrIDs[index] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[index]*WID3;
pos[0] = crd + i;
pos[1] = 2*(j%2) + 1;
pos[2] = 2*(k%2) + 1;
array[vblock::padIndex<PAD>(i_trgt+i,j+PAD,k+PAD)] = vblock::interp_yz<vblock::interpmethod::CIC>(pos,ptr);
}
}
} else { // Neighbor does not exist, return zero values
for (uint32_t i=0; i<PAD; ++i) for (uint32_t k=0; k<WID; ++k) for (uint32_t j=0; j<WID; ++j) {
array[vblock::padIndex<PAD>(i_trgt+i,j+PAD,k+PAD)] = 0.0;
}
}
}
// Copy values from y face neighbors:
for (int j_nbr_off=-1; j_nbr_off<2; j_nbr_off+=2) {
vmesh.getNeighborsExistingAtOffset(blockGID,+0,j_nbr_off,+0,nbrIDs,refLevelDiff);
Real pos[3];
if (j_nbr_off < 0) crd = WID-0.5-(PAD-1);
else crd = 0.5;
uint32_t j_trgt = 0;
if (j_nbr_off > 0) j_trgt = WID+PAD;
if (nbrIDs.size() > 0) {
if (refLevelDiff == -1) { // nbr one level coarser, interpolate
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[0]*WID3;
for (uint32_t j=0; j<PAD; ++j) for (uint32_t k=0; k<WID; ++k) for (uint32_t i=0; i<WID; ++i) {
pos[0] = 2*(i_block%2) + i/2 + 0.5;
pos[1] = crd + j;
pos[2] = 2*(k_block%2) + k/2 + 0.5;
array[vblock::padIndex<PAD>(i+PAD,j_trgt+j,k+PAD)] = vblock::interp_xz<vblock::interpmethod::NGP>(pos,ptr);
}
} else if (refLevelDiff == 0) { // nbr at same level, simple data copy
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[0]*WID3;
uint32_t j_src = 0;
if (j_nbr_off < 0) j_src = WID-PAD;
for (uint32_t j=0; j<PAD; ++j) for (uint32_t k=0; k<WID; ++k) for (uint32_t i=0; i<WID; ++i) {
array[vblock::padIndex<PAD>(i+PAD,j_trgt+j,k+PAD)] = ptr[vblock::index(i,j_src+j,k)];
}
} else if (refLevelDiff == +1) { // nbr one level more refined, interpolate from four neighbors
for (uint32_t j=0; j<PAD; ++j) for (uint32_t k=0; k<WID; ++k) for (uint32_t i=0; i<WID; ++i) {
int index = (k/2)*2 + i/2;
if (nbrIDs[index] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[index]*WID3;
pos[0] = 2*(i%2) + 1;
pos[1] = crd + j;
pos[2] = 2*(k%2) + 1;
array[vblock::padIndex<PAD>(i+PAD,j_trgt+j,k+PAD)] = vblock::interp_xz<vblock::interpmethod::CIC>(pos,ptr);
}
}
} else { // Neighbor does not exist, return zero values
for (uint32_t j=0; j<PAD; ++j) for (uint32_t k=0; k<WID; ++k) for (uint32_t i=0; i<WID; ++i) {
array[vblock::padIndex<PAD>(i+PAD,j_trgt+j,k+PAD)] = 0.0;
}
}
}
// Copy values from z face neighbors:
for (int k_nbr_off=-1; k_nbr_off<2; k_nbr_off+=2) {
vmesh.getNeighborsExistingAtOffset(blockGID,+0,+0,k_nbr_off,nbrIDs,refLevelDiff);
Real pos[3];
uint32_t k_trgt = 0;
if (k_nbr_off > 0) k_trgt = WID+PAD;
if (k_nbr_off < 0) crd = WID-0.5-(PAD-1);
else crd = 0.5;
if (nbrIDs.size() > 0) {
if (refLevelDiff == -1) { // nbr one level coarser, interpolate
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[0]*WID3;
for (uint32_t k=0; k<PAD; ++k) for (uint32_t j=0; j<WID; ++j) for (uint32_t i=0; i<WID; ++i) {
pos[0] = 2*(i_block%2) + i/2 + 0.5;
pos[1] = 2*(j_block%2) + j/2 + 0.5;
pos[2] = crd + k;
array[vblock::padIndex<PAD>(i+PAD,j+PAD,k_trgt+k)] = vblock::interp_xy<vblock::interpmethod::NGP>(pos,ptr);
}
} else if (refLevelDiff == 0) { // nbr at same level, simple data copy
if (nbrIDs[0] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[0]*WID3;
uint32_t k_src = 0;
if (k_nbr_off < 0) k_src = WID-PAD;
for (uint32_t k=0; k<PAD; ++k) for (uint32_t j=0; j<WID; ++j) for (uint32_t i=0; i<WID; ++i) {
array[vblock::padIndex<PAD>(i+PAD,j+PAD,k_trgt+k)] = ptr[vblock::index(i,j,k_src+k)];
}
} else if (refLevelDiff == +1) { // nbr one level more refined, interpolate from four neighbors
for (uint32_t k=0; k<PAD; ++k) for (uint32_t j=0; j<WID; ++j) for (uint32_t i=0; i<WID; ++i) {
int index = (j/2)*2 + i/2;
if (nbrIDs[index] == invalid_local_id()) ptr = null_block_data;
else ptr = src + nbrIDs[index]*WID3;
pos[0] = 2*(i%2) + 1;
pos[1] = 2*(j%2) + 1;
pos[2] = crd + k;
array[vblock::padIndex<PAD>(i+PAD,j+PAD,k_trgt+k)] = vblock::interp_xy<vblock::interpmethod::CIC>(pos,ptr);
}
}
} else { // Neighbor does not exist, return zero values
for (uint32_t k=0; k<PAD; ++k) for (uint32_t j=0; j<WID; ++j) for (uint32_t i=0; i<WID; ++i) {
array[vblock::padIndex<PAD>(i+PAD,j+PAD,k_trgt+k)] = 0.0;
}
}
}
}
inline vmesh::GlobalID SpatialCell::find_velocity_block(uint8_t& refLevel,vmesh::GlobalID cellIndices[3],const int& popID) {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].vmesh.findBlock(refLevel,cellIndices);
}
inline Realf* SpatialCell::get_data(const int& popID) {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].blockContainer.getData();
}
inline const Realf* SpatialCell::get_data(const int& popID) const {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].blockContainer.getData();
}
inline Realf* SpatialCell::get_data(const vmesh::LocalID& blockLID,const int& popID) {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
if (blockLID >= populations[popID].blockContainer.size()) {
std::cerr << "ERROR, block LID out of bounds, blockContainer.size() " << populations[popID].blockContainer.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
if (blockLID == vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>::invalidLocalID()) return null_block_data;
return populations[popID].blockContainer.getData(blockLID);
}
inline const Realf* SpatialCell::get_data(const vmesh::LocalID& blockLID,const int& popID) const {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
if (blockLID >= populations[popID].blockContainer.size()) {
std::cerr << "ERROR, block LID out of bounds, blockContainer.size() " << populations[popID].blockContainer.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
if (blockLID == vmesh::VelocityMesh<vmesh::GlobalID,vmesh::LocalID>::invalidLocalID()) return null_block_data;
return populations[popID].blockContainer.getData(blockLID);
}
inline Real* SpatialCell::get_block_parameters(const int& popID) {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].blockContainer.getParameters();
}
inline const Real* SpatialCell::get_block_parameters(const int& popID) const {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].blockContainer.getParameters();
}
inline Real* SpatialCell::get_block_parameters(const vmesh::LocalID& blockLID,const int& popID) {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
if (blockLID >= populations[popID].blockContainer.size()) {
std::cerr << "ERROR, block LID out of bounds, blockContainer.size() " << populations[popID].blockContainer.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].blockContainer.getParameters(blockLID);
}
inline const Real* SpatialCell::get_block_parameters(const vmesh::LocalID& blockLID,const int& popID) const {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
if (blockLID >= populations[popID].blockContainer.size()) {
std::cerr << "ERROR, block LID out of bounds, blockContainer.size() " << populations[popID].blockContainer.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].blockContainer.getParameters(blockLID);
}
inline Real* SpatialCell::get_cell_parameters() {
return parameters;
}
inline const Real* SpatialCell::get_cell_parameters() const {
return parameters;
}
inline uint8_t SpatialCell::get_maximum_refinement_level(const int& popID) {
return populations[popID].vmesh.getMaxAllowedRefinementLevel();
}
inline vmesh::LocalID SpatialCell::get_number_of_velocity_blocks(const int& popID) const {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].blockContainer.size();
}
/** Get the total number of velocity blocks in this cell, summed over
* all existing particle populations.
* @return Total number of velocity blocks in the cell.*/
inline vmesh::LocalID SpatialCell::get_number_of_all_velocity_blocks() const {
vmesh::LocalID N_blocks = 0;
for (size_t p=0; p<populations.size(); ++p)
N_blocks += populations[p].blockContainer.size();
return N_blocks;
}
inline int SpatialCell::get_number_of_populations() const {
return populations.size();
}
inline const vmesh::LocalID* SpatialCell::get_velocity_grid_length(const int& popID,const uint8_t& refLevel) {
return populations[popID].vmesh.getGridLength(refLevel);
}
inline const Real* SpatialCell::get_velocity_grid_block_size(const int& popID,const uint8_t& refLevel) {
return populations[popID].vmesh.getBlockSize(refLevel);
}
inline const Real* SpatialCell::get_velocity_grid_cell_size(const int& popID,const uint8_t& refLevel) {
return populations[popID].vmesh.getCellSize(refLevel);
}
inline void SpatialCell::get_velocity_block_coordinates(const int& popID,const vmesh::GlobalID& globalID,Real* coords) {
populations[popID].vmesh.getBlockCoordinates(globalID,coords);
}
/*!
Returns the indices of given velocity block
*/
inline velocity_block_indices_t SpatialCell::get_velocity_block_indices(const int& popID,const vmesh::GlobalID block) {
velocity_block_indices_t indices;
uint8_t refLevel;
populations[popID].vmesh.getIndices(block,refLevel,indices[0],indices[1],indices[2]);
return indices;
}
inline velocity_block_indices_t SpatialCell::get_velocity_block_indices(const int& popID,const vmesh::GlobalID block,uint8_t& refLevel) {
velocity_block_indices_t indices;
populations[popID].vmesh.getIndices(block,refLevel,indices[0],indices[1],indices[2]);
return indices;
}
/*!
Returns the velocity block at given indices or error_velocity_block
*/
inline vmesh::GlobalID SpatialCell::get_velocity_block(const int& popID,const velocity_block_indices_t indices,const uint8_t& refLevel) const {
return populations[popID].vmesh.getGlobalID(refLevel,indices[0],indices[1],indices[2]);
}
inline vmesh::GlobalID SpatialCell::get_velocity_block(const int& popID,vmesh::GlobalID blockIndices[3],const uint8_t& refLevel) const {
return populations[popID].vmesh.getGlobalID(refLevel,blockIndices[0],blockIndices[1],blockIndices[2]);
}
/*!
Returns the velocity block at given location or
error_velocity_block if outside of the velocity grid
*/
inline vmesh::GlobalID SpatialCell::get_velocity_block(const int& popID,const Real vx,const Real vy,const Real vz,const uint8_t& refLevel) const {
Real coords[3] = {vx,vy,vz};
return populations[popID].vmesh.getGlobalID(refLevel,coords);
}
inline vmesh::GlobalID SpatialCell::get_velocity_block(const int& popID,const Real* coords,const uint8_t& refLevel) const {
return populations[popID].vmesh.getGlobalID(refLevel,coords);
}
inline vmesh::GlobalID SpatialCell::get_velocity_block_child(const int& popID,const vmesh::GlobalID& blockGID,const uint8_t& refLevel,
const int& i_cell,const int& j_cell,const int& k_cell) {
uint8_t ref = refLevel;
vmesh::LocalID i_child,j_child,k_child;
i_child = 2*i_child + i_cell/2;
j_child = 2*j_child + j_cell/2;
k_child = 2*k_child + k_cell/2;
while (ref != populations[popID].vmesh.getMaxAllowedRefinementLevel()) {
vmesh::LocalID i_child,j_child,k_child;
populations[popID].vmesh.getIndices(blockGID,ref,i_child,j_child,k_child);
return populations[popID].vmesh.getGlobalID(refLevel+1,i_child,j_child,k_child);
}
return populations[popID].vmesh.invalidGlobalID();
}
inline void SpatialCell::get_velocity_block_children_local_ids(
const vmesh::GlobalID& blockGID,
std::vector<vmesh::LocalID>& childrenLIDs,
const int& popID) {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
std::vector<vmesh::GlobalID> childrenGIDs;
populations[popID].vmesh.getChildren(blockGID,childrenGIDs);
childrenLIDs.resize(childrenGIDs.size());
for (size_t c=0; c<childrenGIDs.size(); ++c)
childrenLIDs[c] = populations[popID].vmesh.getLocalID(childrenGIDs[c]);
}
inline vmesh::GlobalID SpatialCell::get_velocity_block_parent(const int& popID,const vmesh::GlobalID& blockGID) {
return populations[popID].vmesh.getParent(blockGID);
}
inline vmesh::GlobalID SpatialCell::get_velocity_block_global_id(const vmesh::LocalID& blockLID,const int& popID) const {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {
std::cerr << "ERROR, popID " << popID << " exceeds populations.size() " << populations.size() << " in ";
std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
exit(1);
}
#endif
return populations[popID].vmesh.getGlobalID(blockLID);
}
inline vmesh::LocalID SpatialCell::get_velocity_block_local_id(const vmesh::GlobalID& blockGID,const int& popID) const {
#ifdef DEBUG_SPATIAL_CELL
if (popID >= populations.size()) {