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
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
// *****************************************************************************
/*!
  \file      src/Inciter/Discretization.cpp
  \copyright 2012-2015 J. Bakosi,
             2016-2018 Los Alamos National Security, LLC.,
             2019-2021 Triad National Security, LLC.
             All rights reserved. See the LICENSE file for details.
  \details   Data and functionality common to all discretization schemes
  \see       Discretization.h and Discretization.C for more info.
*/
// *****************************************************************************

#include "Tags.hpp"
#include "Reorder.hpp"
#include "Vector.hpp"
#include "DerivedData.hpp"
#include "Discretization.hpp"
#include "MeshWriter.hpp"
#include "DiagWriter.hpp"
#include "Inciter/InputDeck/InputDeck.hpp"
#include "Inciter/Options/Scheme.hpp"
#include "Print.hpp"
#include "Around.hpp"
#include "QuinoaBuildConfig.hpp"
#include "ConjugateGradients.hpp"
#include "ALE.hpp"

#ifdef HAS_EXAM2M
  #include "Controller.hpp"
#endif

namespace inciter {

static CkReduction::reducerType PDFMerger;
extern ctr::InputDeck g_inputdeck;
extern ctr::InputDeck g_inputdeck_defaults;

} // inciter::

using inciter::Discretization;

Discretization::Discretization(
  std::size_t meshid,
  const std::vector< CProxy_Discretization >& disc,
  const CProxy_DistFCT& fctproxy,
  const CProxy_ALE& aleproxy,
  const tk::CProxy_ConjugateGradients& conjugategradientsproxy,
  const CProxy_Transporter& transporter,
  const tk::CProxy_MeshWriter& meshwriter,
  const tk::UnsMesh::CoordMap& coordmap,
  const tk::UnsMesh::Chunk& el,
  const tk::CommMaps& msum,
  const std::map< int, std::vector< std::size_t > >& bface,
  const std::vector< std::size_t >& triinpoel,
  int nc ) :
  m_meshid( meshid ),
  m_transfer_complete(),
  m_transfer( g_inputdeck.get< tag::couple, tag::transfer >() ),
  m_disc( disc ),
  m_nchare( nc ),
  m_it( 0 ),
  m_itr( 0 ),
  m_itf( 0 ),
  m_initial( 1 ),
  m_t( g_inputdeck.get< tag::discr, tag::t0 >() ),
  m_lastDumpTime( -std::numeric_limits< tk::real >::max() ),
  m_physFieldFloor( 0.0 ),
  m_physHistFloor( 0.0 ),
  m_rangeFieldFloor(
    g_inputdeck.get< tag::output, tag::range, tag::field >().size(), 0.0 ),
  m_rangeHistFloor(
    g_inputdeck.get< tag::output, tag::range, tag::history >().size(), 0.0 ),
  m_dt( g_inputdeck.get< tag::discr, tag::dt >() ),
  m_dtn( m_dt ),
  m_nvol( 0 ),
  m_fct( fctproxy ),
  m_ale( aleproxy ),
  m_transporter( transporter ),
  m_meshwriter( meshwriter ),
  m_el( el ),     // fills m_inpoel, m_gid, m_lid
  m_coord( setCoord( coordmap ) ),
  m_coordn( m_coord ),
  m_nodeCommMap(),
  m_edgeCommMap(),
  m_meshvol( 0.0 ),
  m_v( m_gid.size(), 0.0 ),
  m_vol( m_gid.size(), 0.0 ),
  m_volc(),
  m_voln( m_vol ),
  m_vol0( m_inpoel.size()/4, 0.0 ),
  m_bid(),
  m_timer(),
  m_refined( 0 ),
  m_prevstatus( std::chrono::high_resolution_clock::now() ),
  m_nrestart( 0 ),
  m_histdata(),
  m_nsrc( 0 ),
  m_ndst( 0 ),
  m_meshvel( 0, 3 ),
  m_meshvel_converged( true ),
  m_bface( bface ),
  m_triinpoel( triinpoel )
// *****************************************************************************
//  Constructor
//! \param[in] meshid Mesh ID
//! \param[in] disc All Discretization proxies (one per mesh)
//! \param[in] fctproxy Distributed FCT proxy
//! \param[in] aleproxy Distributed ALE proxy
//! \param[in] conjugategradientsproxy Distributed Conjugrate Gradients linear
//!   solver proxy
//! \param[in] transporter Host (Transporter) proxy
//! \param[in] meshwriter Mesh writer proxy
//! \param[in] coordmap Coordinates of mesh nodes and their global IDs
//! \param[in] msum Communication maps associated to chare IDs bordering the
//!   mesh chunk we operate on
//! \param[in] nc Total number of Discretization chares
// *****************************************************************************
{
  Assert( !m_inpoel.empty(), "No elements assigned to Discretization chare" );
  Assert( tk::positiveJacobians( m_inpoel, m_coord ),
          "Jacobian in input mesh to Discretization non-positive" );
  #if not defined(__INTEL_COMPILER) || defined(NDEBUG)
  // The above ifdef skips running the conformity test with the intel compiler
  // in debug mode only. This is necessary because in tk::conforming(), filling
  // up the map can fail with some meshes (only in parallel), e.g., tube.exo,
  // used by some regression tests, due to the intel compiler generating some
  // garbage incorrect code - only in debug, only in parallel, only with that
  // mesh.
  Assert( tk::conforming( m_inpoel, m_coord ),
          "Input mesh to Discretization not conforming" );
  #endif

  // Store communication maps
  for (const auto& [ c, maps ] : msum) {
    m_nodeCommMap[c] = maps.get< tag::node >();
    m_edgeCommMap[c] = maps.get< tag::edge >();
  }

  // Get ready for computing/communicating nodal volumes
  startvol();

  // Count the number of mesh nodes at which we receive data from other chares
  // and compute map associating boundary-chare node ID to global node ID
  std::vector< std::size_t > c( tk::sumvalsize( m_nodeCommMap ) );
  std::size_t j = 0;
  for (const auto& [ch,n] : m_nodeCommMap) for (auto i : n) c[j++] = i;
  tk::unique( c );
  m_bid = tk::assignLid( c );

  // Find host elements of user-specified points where time histories are
  // saved, and save the shape functions evaluated at the point locations
  const auto& pt = g_inputdeck.get< tag::history, tag::point >();
  const auto& id = g_inputdeck.get< tag::history, tag::id >();
  for (std::size_t p=0; p<pt.size(); ++p) {
    std::array< tk::real, 4 > N;
    const auto& l = pt[p];
    for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
      if (tk::intet( m_coord, m_inpoel, l, e, N )) {
        m_histdata.push_back( HistData{{ id[p], e, {l[0],l[1],l[2]}, N }} );
        break;
      }
    }
  }

  // Insert DistFCT chare array element if FCT is needed. Note that even if FCT
  // is configured false in the input deck, at this point, we still need the FCT
  // object as FCT is still being performed, only its results are ignored.
  const auto sch = g_inputdeck.get< tag::discr, tag::scheme >();
  const auto nprop = g_inputdeck.get< tag::component >().nprop();
  if (sch == ctr::SchemeType::DiagCG)
    m_fct[ thisIndex ].insert( m_nchare, m_gid.size(), nprop,
                               m_nodeCommMap, m_bid, m_lid, m_inpoel );

  // Insert ConjugrateGradients solver chare array element if needed
  if (g_inputdeck.get< tag::ale, tag::ale >()) {
    m_ale[ thisIndex ].insert( conjugategradientsproxy,
                               m_coord, m_inpoel,
                               m_gid, m_lid, m_nodeCommMap );
  } else {
    m_meshvel.resize( m_gid.size() );
  }

  // Register mesh with mesh-transfer lib
  if (m_disc.size() == 1 || m_transfer.empty()) {
    // skip transfer if single mesh or if not involved in coupling
    transferInit();
  } else {
    #ifdef HAS_EXAM2M
    if (thisIndex == 0) {
      exam2m::addMesh( thisProxy, m_nchare,
        CkCallback( CkIndex_Discretization::transferInit(), thisProxy ) );
      std::cout << "Disc: " << m_meshid << " m2m::addMesh()\n";
    }
    #else
    transferInit();
    #endif
  }
}

void
Discretization::transferInit()
// *****************************************************************************
// Our mesh has been registered with the mesh-to-mesh transfer library (if
// coupled to other solver)
// *****************************************************************************
{
  // Compute number of mesh points owned
  std::size_t npoin = m_gid.size();
  for (auto g : m_gid) if (tk::slave(m_nodeCommMap,g,thisIndex)) --npoin;

  // Tell the RTS that the Discretization chares have been created and compute
  // the total number of mesh points across the distributed mesh
  std::vector< std::size_t > meshdata{ m_meshid, npoin };
  contribute( meshdata, CkReduction::sum_ulong,
    CkCallback( CkReductionTarget(Transporter,disccreated), m_transporter ) );
}

void
Discretization::meshvelStart(
  const tk::UnsMesh::Coords vel,
  const std::vector< tk::real >& soundspeed,
  const std::unordered_map< int,
    std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& bnorm,
  tk::real adt,
  CkCallback done ) const
// *****************************************************************************
// Start computing new mesh velocity for ALE mesh motion
//! \param[in] vel Fluid velocity at mesh nodes
//! \param[in] soundspeed Speed of sound at mesh nodes
//! \param[in] bnorm Face normals in boundary points associated to side sets
//! \param[in] adt alpha*dt of the RK time step
//! \param[in] done Function to continue with when mesh velocity has been
//!   computed
// *****************************************************************************
{
  if (g_inputdeck.get< tag::ale, tag::ale >())
    m_ale[ thisIndex ].ckLocal()->start( vel, soundspeed, done,
      m_coord, m_coordn, m_vol0, m_vol, bnorm, m_initial, m_it, m_t, adt );
  else
    done.send();
}

const tk::Fields&
Discretization::meshvel() const
// *****************************************************************************
//! Query the mesh velocity
//! \return Mesh velocity
// *****************************************************************************
{
  if (g_inputdeck.get< tag::ale, tag::ale >())
    return m_ale[ thisIndex ].ckLocal()->meshvel();
  else
    return m_meshvel;
}

void
Discretization::meshvelBnd(
  const std::map< int, std::vector< std::size_t > >& bface,
  const std::map< int, std::vector< std::size_t > >& bnode,
  const std::vector< std::size_t >& triinpoel) const
// *****************************************************************************
// Query ALE mesh velocity boundary condition node lists and node lists at
// which ALE moves boundaries
// *****************************************************************************
{
  if (g_inputdeck.get< tag::ale, tag::ale >())
    m_ale[ thisIndex ].ckLocal()->meshvelBnd( bface, bnode, triinpoel );
}

void
Discretization::meshvelConv()
// *****************************************************************************
//! Assess and record mesh velocity linear solver convergence
// *****************************************************************************
{
  auto smoother = g_inputdeck.get< tag::ale, tag::smoother >();

  if (g_inputdeck.get< tag::ale, tag::ale >() &&
      (smoother == ctr::MeshVelocitySmootherType::LAPLACE or
       smoother == ctr::MeshVelocitySmootherType::HELMHOLTZ))
  {
    m_meshvel_converged &= m_ale[ thisIndex ].ckLocal()->converged();
  }
}

void
Discretization::transferCallback( std::vector< CkCallback >& cb )
// *****************************************************************************
// Receive a list of callbacks from our own child solver
//! \param[in] cb List of callbacks
//! \details This is called by our child solver, either when it is coupled to
//!    another solver or not.
// *****************************************************************************
{
  // Store callback for when there is no transfer we are involved in
  m_transfer_complete = cb.back();
  cb.pop_back();

  // Distribute callbacks
  for (auto& t : m_transfer) {
    // If we are a source of a transfer, send callback to the destination solver
    if (m_meshid == t.src) {
      Assert( !cb.empty(), "Insufficient number of src callbacks, meshid: " +
                           std::to_string(m_meshid) );
      m_disc[ t.dst ][ thisIndex ].comcb( m_meshid, cb.back() );
      cb.pop_back();
    // If we are a destination of a callback, store it
    } else if (m_meshid == t.dst) {
      Assert( !cb.empty(), "Insufficient number of dst callbacks, meshid: " +
                           std::to_string(m_meshid) );
      t.cb.push_back( cb.back() );
      cb.pop_back();
      //t.cbs.push_back( m_meshid );    // only for debugging
    }
  }
  Assert( cb.empty(), "Not all callbacks have been processed" );

  if (transferCallbacksComplete()) comfinal();
}

void
Discretization::comcb( std::size_t srcmeshid, CkCallback c )
// *****************************************************************************
// Receive mesh transfer callbacks from source mesh/solver
//! \param[in] srcmeshid Source mesh (solver) id
//! \param[in] c Callback received
// *****************************************************************************
{
  // Store received mesh transfer callback from source mesh/solver
  for (auto& t : m_transfer)
    if (srcmeshid == t.src && m_meshid == t.dst) {
      t.cb.push_back( c );
      //t.cbs.push_back( srcmeshid );   // only for debugging
    }

  if (transferCallbacksComplete()) comfinal();
}

bool
Discretization::transferCallbacksComplete() const
// *****************************************************************************
// Determine if communication of mesh transfer callbacks is complete
//! \return True if communication of mesh transfer callbacks have been
//!   completed on this solver
// *****************************************************************************
{
  bool c = true;

  // Our callbacks are complete if all transfers we are involved in as a
  // destination have exactly two callbacks.
  for (const auto& t : m_transfer)
    if (m_meshid == t.dst && t.cb.size() != 2)
      c = false;<--- Consider using std::any_of, std::all_of, std::none_of, or std::accumulate algorithm instead of a raw loop.

  return c;
}

void
Discretization::comfinal()
// *****************************************************************************
// Finish setting up communication maps and solution transfer callbacks
// *****************************************************************************
{
//  std::cout << "m:" << m_meshid << ": transfer: ";
//  for (const auto& t : m_transfer) {
//    std::cout << t.src << "->" << t.dst << ' ';
//    if (t.cb.size() > 0) {
//      std::cout << "cb: ";
//      for (auto m : t.cbs) std::cout << m << ' ';
//    }
//  }
//  std::cout << '\n';

  // Generate own subset of solver/mesh transfer list
  for (const auto& t : m_transfer)
    if (t.src == m_meshid || t.dst == m_meshid)
      m_mytransfer.push_back( t );<--- Consider using std::copy_if algorithm instead of a raw loop.

//  std::cout << "m:" << m_meshid << ": mytransfer: ";
//  for (const auto& t : m_mytransfer) {
//    std::cout << t.src << "->" << t.dst << ' ';
//    if (t.cb.size() > 0) {
//      std::cout << "cb: ";
//      for (auto m : t.cbs) std::cout << m << ' ';
//    }
//  }
//  std::cout << '\n';

  // Signal the runtime system that the workers have been created
  std::vector< std::size_t > meshdata{ /* initial */ 1, m_meshid };
  contribute( meshdata, CkReduction::sum_ulong,
    CkCallback(CkReductionTarget(Transporter,comfinal), m_transporter) );
}

void
Discretization::transfer( [[maybe_unused]] const tk::Fields& u )
// *****************************************************************************
//  Start solution transfer (if coupled)
//! \param[in] u Solution to transfer from/to
// *****************************************************************************
{
  if (m_mytransfer.empty()) {   // skip transfer if not involved in coupling
    m_transfer_complete.send();
  } else {
    // Pass source and destination meshes to mesh transfer lib (if coupled)
    #ifdef HAS_EXAM2M
    Assert( m_nsrc < m_mytransfer.size(), "Indexing out of mytransfer[src]" );
    if (m_mytransfer[m_nsrc].src == m_meshid) {
      exam2m::setSourceTets( thisProxy, thisIndex, &m_inpoel, &m_coord, u );
      ++m_nsrc;
      //std::cout << m_meshid << " src\n";
    } else {
      m_nsrc = 0;
    }
    Assert( m_ndst < m_mytransfer.size(), "Indexing out of mytransfer[dst]" );
    if (m_mytransfer[m_ndst].dst == m_meshid) {
      exam2m::setDestPoints( thisProxy, thisIndex, &m_coord, u,
                             m_mytransfer[m_ndst].cb );
      ++m_ndst;
      //std::cout << m_meshid << " dst\n";
    } else {
      m_ndst = 0;
    }
    #else
    m_transfer_complete.send();
    #endif
  }
}

std::vector< std::size_t >
Discretization::bndel() const
// *****************************************************************************
// Find elements along our mesh chunk boundary
//! \return List of local element ids that have at least a single node
//!   contributing to a chare boundary
// *****************************************************************************
{
  // Lambda to find out if a mesh node is shared with another chare
  auto shared = [this]( std::size_t i ){
    for (const auto& [c,n] : m_nodeCommMap)
      if (n.find(i) != end(n)) return true;
    return false;
  };

  // Find elements along our mesh chunk boundary
  std::vector< std::size_t > e;
  for (std::size_t n=0; n<m_inpoel.size(); ++n)
    if (shared( m_gid[ m_inpoel[n] ] )) e.push_back( n/4 );
  tk::unique( e );

  return e;
}

void
Discretization::resizePostAMR( const tk::UnsMesh::Chunk& chunk,
                               const tk::UnsMesh::Coords& coord,
                               const tk::NodeCommMap& nodeCommMap )
// *****************************************************************************
//  Resize mesh data structures after mesh refinement
//! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
//! \param[in] coord New mesh node coordinates
//! \param[in] nodeCommMap New node communication map
// *****************************************************************************
{
  m_el = chunk;         // updates m_inpoel, m_gid, m_lid
  m_nodeCommMap = nodeCommMap;        // update node communication map

  // Update mesh volume container size
  m_vol.resize( m_gid.size(), 0.0 );

  // Generate local ids for new chare boundary global ids
  std::size_t bid = m_bid.size();
  for (const auto& [ neighborchare, sharednodes ] : m_nodeCommMap)
    for (auto g : sharednodes)
      if (m_bid.find(g) == end(m_bid))
        m_bid[g] = bid++;<--- Searching before insertion is not necessary.

  // update mesh node coordinates
  m_coord = coord;

  // we are no longer during setup
  m_initial = 0;
}

void
Discretization::startvol()
// *****************************************************************************
//  Get ready for (re-)computing/communicating nodal volumes
// *****************************************************************************
{
  m_nvol = 0;
  thisProxy[ thisIndex ].wait4vol();

  // Zero out mesh volume container
  std::fill( begin(m_vol), end(m_vol), 0.0 );

  // Clear receive buffer that will be used for collecting nodal volumes
  m_volc.clear();
}

void
Discretization::registerReducers()
// *****************************************************************************
//  Configure Charm++ reduction types
//!  \details Since this is a [initnode] routine, see the .ci file, the
//!   Charm++ runtime system executes the routine exactly once on every
//!   logical node early on in the Charm++ init sequence. Must be static as
//!   it is called without an object. See also: Section "Initializations at
//!   Program Startup" at in the Charm++ manual
//!   http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
// *****************************************************************************
{
  PDFMerger = CkReduction::addReducer( tk::mergeUniPDFs );
}

tk::UnsMesh::Coords
Discretization::setCoord( const tk::UnsMesh::CoordMap& coordmap )
// *****************************************************************************
// Set mesh coordinates based on coordinates map
// *****************************************************************************
{
  Assert( coordmap.size() == m_gid.size(), "Size mismatch" );
  Assert( coordmap.size() == m_lid.size(), "Size mismatch" );

  tk::UnsMesh::Coords coord;
  coord[0].resize( coordmap.size() );
  coord[1].resize( coordmap.size() );
  coord[2].resize( coordmap.size() );

  for (const auto& [ gid, coords ] : coordmap) {
    auto i = tk::cref_find( m_lid, gid );
    coord[0][i] = coords[0];
    coord[1][i] = coords[1];
    coord[2][i] = coords[2];
  }

  return coord;
}

void
Discretization::remap(
  const std::unordered_map< std::size_t, std::size_t >& map )
// *****************************************************************************
//  Remap mesh data based on new local ids
//! \param[in] map Mapping of old->new local ids
// *****************************************************************************
{
  // Remap connectivity containing local IDs
  for (auto& l : m_inpoel) l = tk::cref_find(map,l);

  // Remap global->local id map
  for (auto& [g,l] : m_lid) l = tk::cref_find(map,l);

  // Remap global->local id map
  auto maxid = std::numeric_limits< std::size_t >::max();
  std::vector< std::size_t > newgid( m_gid.size(), maxid );
  for (const auto& [o,n] : map) newgid[n] = m_gid[o];
  m_gid = std::move( newgid );

  Assert( std::all_of( m_gid.cbegin(), m_gid.cend(),
            [=](std::size_t i){ return i < maxid; } ),
          "Not all gid have been remapped" );

  // Remap nodal volumes (with contributions along chare-boundaries)
  std::vector< tk::real > newvol( m_vol.size(), 0.0 );
  for (const auto& [o,n] : map) newvol[n] = m_vol[o];
  m_vol = std::move( newvol );

  // Remap nodal volumes (without contributions along chare-boundaries)
  std::vector< tk::real > newv( m_v.size(), 0.0 );
  for (const auto& [o,n] : map) newv[n] = m_v[o];
  m_v = std::move( newv );

  // Remap locations of node coordinates
  tk::UnsMesh::Coords newcoord;
  auto npoin = m_coord[0].size();
  newcoord[0].resize( npoin );
  newcoord[1].resize( npoin );
  newcoord[2].resize( npoin );
  for (const auto& [o,n] : map) {
    newcoord[0][n] = m_coord[0][o];
    newcoord[1][n] = m_coord[1][o];
    newcoord[2][n] = m_coord[2][o];
  }
  m_coord = std::move( newcoord );
}

void
Discretization::setRefiner( const CProxy_Refiner& ref )
// *****************************************************************************
//  Set Refiner Charm++ proxy
//! \param[in] ref Incoming refiner proxy to store
// *****************************************************************************
{
  m_refiner = ref;
}

void
Discretization::vol()
// *****************************************************************************
// Sum mesh volumes to nodes, start communicating them on chare-boundaries
// *****************************************************************************
{
  const auto& x = m_coord[0];
  const auto& y = m_coord[1];
  const auto& z = m_coord[2];

  // Compute nodal volumes on our chunk of the mesh
  for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
    const std::array< std::size_t, 4 > N{{ m_inpoel[e*4+0], m_inpoel[e*4+1],
                                           m_inpoel[e*4+2], m_inpoel[e*4+3] }};
    // compute element Jacobi determinant * 5/120 = element volume / 4
    const std::array< tk::real, 3 >
      ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
      ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
      da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
    const auto J = tk::triple( ba, ca, da ) * 5.0 / 120.0;
    ErrChk( J > 0, "Element Jacobian non-positive: PE:" +
                   std::to_string(CkMyPe()) + ", node IDs: " +
                   std::to_string(m_gid[N[0]]) + ',' +
                   std::to_string(m_gid[N[1]]) + ',' +
                   std::to_string(m_gid[N[2]]) + ',' +
                   std::to_string(m_gid[N[3]]) + ", coords: (" +
                   std::to_string(x[N[0]]) + ", " +
                   std::to_string(y[N[0]]) + ", " +
                   std::to_string(z[N[0]]) + "), (" +
                   std::to_string(x[N[1]]) + ", " +
                   std::to_string(y[N[1]]) + ", " +
                   std::to_string(z[N[1]]) + "), (" +
                   std::to_string(x[N[2]]) + ", " +
                   std::to_string(y[N[2]]) + ", " +
                   std::to_string(z[N[2]]) + "), (" +
                   std::to_string(x[N[3]]) + ", " +
                   std::to_string(y[N[3]]) + ", " +
                   std::to_string(z[N[3]]) + ')' );
    // scatter add V/4 to nodes
    for (std::size_t j=0; j<4; ++j) m_vol[N[j]] += J;

    // save element volumes at t=t0
    if (m_it == 0) m_vol0[e] = J * 4.0;
  }

  // Store nodal volumes without contributions from other chares on
  // chare-boundaries
  m_v = m_vol;

  // Send our nodal volume contributions to neighbor chares
  if (m_nodeCommMap.empty())
   totalvol();
  else
    for (const auto& [c,n] : m_nodeCommMap) {
      std::vector< tk::real > v( n.size() );
      std::size_t j = 0;
      for (auto i : n) v[ j++ ] = m_vol[ tk::cref_find(m_lid,i) ];
      thisProxy[c].comvol( std::vector<std::size_t>(begin(n), end(n)), v );
    }

  ownvol_complete();
}

void
Discretization::comvol( const std::vector< std::size_t >& gid,
                        const std::vector< tk::real >& nodevol )
// *****************************************************************************
//  Receive nodal volumes on chare-boundaries
//! \param[in] gid Global mesh node IDs at which we receive volume contributions
//! \param[in] nodevol Partial sums of nodal volume contributions to
//!    chare-boundary nodes
//! \details This function receives contributions to m_vol, which stores the
//!   nodal volumes. While m_vol stores own contributions, m_volc collects the
//!   neighbor chare contributions during communication. This way work on m_vol
//!   and m_volc is overlapped. The contributions are applied in totalvol().
// *****************************************************************************
{
  Assert( nodevol.size() == gid.size(), "Size mismatch" );

  for (std::size_t i=0; i<gid.size(); ++i)
    m_volc[ gid[i] ] += nodevol[i];

  if (++m_nvol == m_nodeCommMap.size()) {
    m_nvol = 0;
    comvol_complete();
  }
}

void
Discretization::totalvol()
// *****************************************************************************
// Sum mesh volumes and contribute own mesh volume to total volume
// *****************************************************************************
{
  // Add received contributions to nodal volumes
  for (const auto& [gid, vol] : m_volc)
    m_vol[ tk::cref_find(m_lid,gid) ] += vol;

  // Clear receive buffer
  tk::destroy(m_volc);

  // Sum mesh volume to host
  std::vector< tk::real > tvol{ 0.0,
                                static_cast<tk::real>(m_initial),
                                static_cast<tk::real>(m_meshid) };
  for (auto v : m_v) tvol[0] += v;
  contribute( tvol, CkReduction::sum_double,
    CkCallback(CkReductionTarget(Transporter,totalvol), m_transporter) );
}

void
Discretization::stat( tk::real mesh_volume )
// *****************************************************************************
// Compute mesh cell statistics
//! \param[in] mesh_volume Total mesh volume
// *****************************************************************************
{
  // Store total mesh volume
  m_meshvol = mesh_volume;

  const auto& x = m_coord[0];
  const auto& y = m_coord[1];
  const auto& z = m_coord[2];

  auto MIN = -std::numeric_limits< tk::real >::max();
  auto MAX = std::numeric_limits< tk::real >::max();
  std::vector< tk::real > min{ MAX, MAX, MAX };
  std::vector< tk::real > max{ MIN, MIN, MIN };
  std::vector< tk::real > sum{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  tk::UniPDF edgePDF( 1e-4 );
  tk::UniPDF volPDF( 1e-4 );
  tk::UniPDF ntetPDF( 1e-4 );

  // Compute points surrounding points
  auto psup = tk::genPsup( m_inpoel, 4, tk::genEsup(m_inpoel,4) );
  Assert( psup.second.size()-1 == m_gid.size(),
          "Number of mesh points and number of global IDs unequal" );

  // Compute edge length statistics
  // Note that while the min and max edge lengths are independent of the number
  // of CPUs (by the time they are aggregated across all chares), the sum of
  // the edge lengths and the edge length PDF are not. This is because the
  // edges on the chare-boundary are counted multiple times and we
  // conscientiously do not make an effort to precisely compute this, because
  // that would require communication and more complex logic. Since these
  // statistics are intended as simple average diagnostics, we ignore these
  // small differences. For reproducible average edge lengths and edge length
  // PDFs, run the mesh in serial.
  for (std::size_t p=0; p<m_gid.size(); ++p)
    for (auto i : tk::Around(psup,p)) {
       const auto dx = x[ i ] - x[ p ];
       const auto dy = y[ i ] - y[ p ];
       const auto dz = z[ i ] - z[ p ];
       const auto length = std::sqrt( dx*dx + dy*dy + dz*dz );
       if (length < min[0]) min[0] = length;
       if (length > max[0]) max[0] = length;
       sum[0] += 1.0;
       sum[1] += length;
       edgePDF.add( length );
    }

  // Compute mesh cell volume statistics
  for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
    const std::array< std::size_t, 4 > N{{ m_inpoel[e*4+0], m_inpoel[e*4+1],
                                           m_inpoel[e*4+2], m_inpoel[e*4+3] }};
    const std::array< tk::real, 3 >
      ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
      ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
      da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
    const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
    if (L < min[1]) min[1] = L;
    if (L > max[1]) max[1] = L;
    sum[2] += 1.0;
    sum[3] += L;
    volPDF.add( L );
  }

  // Contribute stats of number of tetrahedra (ntets)
  sum[4] = 1.0;
  min[2] = max[2] = sum[5] = static_cast< tk::real >( m_inpoel.size() / 4 );
  ntetPDF.add( min[2] );

  min.push_back( static_cast<tk::real>(m_meshid) );
  max.push_back( static_cast<tk::real>(m_meshid) );
  sum.push_back( static_cast<tk::real>(m_meshid) );

  // Contribute to mesh statistics across all Discretization chares
  contribute( min, CkReduction::min_double,
    CkCallback(CkReductionTarget(Transporter,minstat), m_transporter) );
  contribute( max, CkReduction::max_double,
    CkCallback(CkReductionTarget(Transporter,maxstat), m_transporter) );
  contribute( sum, CkReduction::sum_double,
    CkCallback(CkReductionTarget(Transporter,sumstat), m_transporter) );

  // Serialize PDFs to raw stream
  auto stream = tk::serialize( m_meshid, { edgePDF, volPDF, ntetPDF } );
  // Create Charm++ callback function for reduction of PDFs with
  // Transporter::pdfstat() as the final target where the results will appear.
  CkCallback cb( CkIndex_Transporter::pdfstat(nullptr), m_transporter );
  // Contribute serialized PDF of partial sums to host via Charm++ reduction
  contribute( stream.first, stream.second.get(), PDFMerger, cb );
}

void
Discretization::boxvol(
  const std::vector< std::unordered_set< std::size_t > >& nodes )
// *****************************************************************************
// Compute total box IC volume
//! \param[in] nodes Node list contributing to box IC volume (for each IC box)
// *****************************************************************************
{
  // Compute partial box IC volume (just add up all boxes)
  tk::real boxvol = 0.0;
  for (const auto& b : nodes) for (auto i : b) boxvol += m_v[i];<--- Consider using std::accumulate algorithm instead of a raw loop.

  // Sum up box IC volume across all chares
  std::vector< tk::real > meshdata{ boxvol, static_cast<tk::real>(m_meshid) };
  contribute( meshdata, CkReduction::sum_double,
    CkCallback(CkReductionTarget(Transporter,boxvol), m_transporter) );
}

void
Discretization::write(
  const std::vector< std::size_t >& inpoel,
  const tk::UnsMesh::Coords& coord,
  const std::map< int, std::vector< std::size_t > >& bface,
  const std::map< int, std::vector< std::size_t > >& bnode,
  const std::vector< std::size_t >& triinpoel,
  const std::vector< std::string>& elemfieldnames,
  const std::vector< std::string>& nodefieldnames,
  const std::vector< std::string>& nodesurfnames,
  const std::vector< std::vector< tk::real > >& elemfields,
  const std::vector< std::vector< tk::real > >& nodefields,
  const std::vector< std::vector< tk::real > >& nodesurfs,
  CkCallback c )
// *****************************************************************************
//  Output mesh and fields data (solution dump) to file(s)
//! \param[in] inpoel Mesh connectivity for the mesh chunk to be written
//! \param[in] coord Node coordinates of the mesh chunk to be written
//! \param[in] bface Map of boundary-face lists mapped to corresponding side set
//!   ids for this mesh chunk
//! \param[in] bnode Map of boundary-node lists mapped to corresponding side set
//!   ids for this mesh chunk
//! \param[in] triinpoel Interconnectivity of points and boundary-face in this
//!   mesh chunk
//! \param[in] elemfieldnames Names of element fields to be output to file
//! \param[in] nodefieldnames Names of node fields to be output to file
//! \param[in] nodesurfnames Names of node surface fields to be output to file
//! \param[in] elemfields Field data in mesh elements to output to file
//! \param[in] nodefields Field data in mesh nodes to output to file
//! \param[in] nodesurfs Surface field data in mesh nodes to output to file
//! \param[in] c Function to continue with after the write
//! \details Since m_meshwriter is a Charm++ chare group, it never migrates and
//!   an instance is guaranteed on every PE. We index the first PE on every
//!   logical compute node. In Charm++'s non-SMP mode, a node is the same as a
//!   PE, so the index is the same as CkMyPe(). In SMP mode the index is the
//!   first PE on every logical node. In non-SMP mode this yields one or more
//!   output files per PE with zero or non-zero virtualization, respectively. If
//!   there are multiple chares on a PE, the writes are serialized per PE, since
//!   only a single entry method call can be executed at any given time. In SMP
//!   mode, still the same number of files are output (one per chare), but the
//!   output is serialized through the first PE of each compute node. In SMP
//!   mode, channeling multiple files via a single PE on each node is required
//!   by NetCDF and HDF5, as well as ExodusII, since none of these libraries are
//!   thread-safe.
// *****************************************************************************
{
  // If the previous iteration refined (or moved) the mesh or this is called
  // before the first time step, we also output the mesh.
  bool meshoutput = m_itf == 0 ? true : false;

  auto eps = std::numeric_limits< tk::real >::epsilon();
  bool fieldoutput = false;

  // Output field data only if there is no dump at this physical time yet
  if (std::abs(m_lastDumpTime - m_t) > eps ) {
    m_lastDumpTime = m_t;
    ++m_itf;
    fieldoutput = true;
  }

  m_meshwriter[ CkNodeFirst( CkMyNode() ) ].
    write( m_meshid, meshoutput, fieldoutput, m_itr, m_itf, m_t, thisIndex,
           g_inputdeck.get< tag::cmd, tag::io, tag::output >(),
           inpoel, coord, bface, bnode, triinpoel, elemfieldnames,
           nodefieldnames, nodesurfnames, elemfields, nodefields, nodesurfs,
           g_inputdeck.outsets(), c );
}

void
Discretization::setdt( tk::real newdt )
// *****************************************************************************
// Set time step size
//! \param[in] newdt Size of the new time step
// *****************************************************************************
{
  m_dtn = m_dt;
  m_dt = newdt;

  // Truncate the size of last time step
  const auto term = g_inputdeck.get< tag::discr, tag::term >();
  if (m_t+m_dt > term) m_dt = term - m_t;
}

void
Discretization::next()
// *****************************************************************************
// Prepare for next step
// *****************************************************************************
{
  // Update floor of physics time divided by output interval times
  const auto eps = std::numeric_limits< tk::real >::epsilon();
  const auto ft = g_inputdeck.get< tag::output, tag::time, tag::field >();
  if (ft > eps) m_physFieldFloor = std::floor( m_t / ft );
  const auto ht = g_inputdeck.get< tag::output, tag::time, tag::history >();
  if (ht > eps) m_physHistFloor = std::floor( m_t / ht );

  // Update floors of physics time divided by output interval times for ranges
  const auto& rf = g_inputdeck.get< tag::output, tag::range, tag::field >();
  for (std::size_t i=0; i<rf.size(); ++i)
    if (m_t > rf[i][0] and m_t < rf[i][1])
      m_rangeFieldFloor[i] = std::floor( m_t / rf[i][2] );
  const auto& rh = g_inputdeck.get< tag::output, tag::range, tag::history >();
  for (std::size_t i=0; i<rh.size(); ++i)
    if (m_t > rh[i][0] and m_t < rh[i][1])
      m_rangeHistFloor[i] = std::floor( m_t / rh[i][2] );

  ++m_it;
  m_t += m_dt;
}

void
Discretization::grindZero()
// *****************************************************************************
//  Zero grind-time
// *****************************************************************************
{
  m_prevstatus = std::chrono::high_resolution_clock::now();

  if (thisIndex == 0 && m_meshid == 0) {
    const auto verbose = g_inputdeck.get< tag::cmd, tag::verbose >();
    const auto& def =
      g_inputdeck_defaults.get< tag::cmd, tag::io, tag::screen >();
    tk::Print print( g_inputdeck.get< tag::cmd >().logname( def, m_nrestart ),
                     verbose ? std::cout : std::clog,
                     std::ios_base::app );
    print.diag( "Starting time stepping ..." );
  }
}

bool
Discretization::restarted( int nrestart )
// *****************************************************************************
//  Detect if just returned from a checkpoint and if so, zero timers
//! \param[in] nrestart Number of times restarted
//! \return True if restart detected
// *****************************************************************************
{
  // Detect if just restarted from checkpoint:
  //   nrestart == -1 if there was no checkpoint this step
  //   d->Nrestart() == nrestart if there was a checkpoint this step
  //   if both false, just restarted from a checkpoint
  bool restarted = nrestart != -1 and m_nrestart != nrestart;

   // If just restarted from checkpoint
  if (restarted) {
    // Update number of restarts
    m_nrestart = nrestart;
    // Start timer measuring time stepping wall clock time
    m_timer.zero();
    // Zero grind-timer
    grindZero();
  }

  return restarted;
}

std::string
Discretization::histfilename( const std::string& id,
                              kw::precision::info::expect::type precision )
// *****************************************************************************
//  Construct history output filename
//! \param[in] id History point id
//! \param[in] precision Floating point precision to use for output
//! \return History file name
// *****************************************************************************
{
  auto of = g_inputdeck.get< tag::cmd, tag::io, tag::output >();
  std::stringstream ss;

  ss << std::setprecision(static_cast<int>(precision)) << of << ".hist." << id;

  return ss.str();
}

void
Discretization::histheader( std::vector< std::string >&& names )
// *****************************************************************************
//  Output headers for time history files (one for each point)
//! \param[in] names History output variable names
// *****************************************************************************
{
  for (const auto& h : m_histdata) {
    auto prec = g_inputdeck.get< tag::prec, tag::history >();
    tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
                       g_inputdeck.get< tag::flformat, tag::history >(),
                       prec );
    hw.header( names );
  }
}

void
Discretization::history( std::vector< std::vector< tk::real > >&& data )
// *****************************************************************************
//  Output time history for a time step
//! \param[in] data Time history data for all variables and equations integrated
// *****************************************************************************
{
  Assert( data.size() == m_histdata.size(), "Size mismatch" );

  std::size_t i = 0;
  for (const auto& h : m_histdata) {
    auto prec = g_inputdeck.get< tag::prec, tag::history >();
    tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
                       g_inputdeck.get< tag::flformat, tag::history >(),
                       prec,
                       std::ios_base::app );
    hw.diag( m_it, m_t, m_dt, data[i] );
    ++i;
  }
}

bool
Discretization::fielditer() const
// *****************************************************************************
//  Decide if field output iteration count interval is hit
//! \return True if field output iteration count interval is hit
// *****************************************************************************
{
  if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;

  return m_it % g_inputdeck.get< tag::output, tag::iter, tag::field >() == 0;
}

bool
Discretization::fieldtime() const
// *****************************************************************************
//  Decide if field output physics time interval is hit
//! \return True if field output physics time interval is hit
// *****************************************************************************
{
  if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;

  const auto eps = std::numeric_limits< tk::real >::epsilon();
  const auto ft = g_inputdeck.get< tag::output, tag::time, tag::field >();

  if (ft < eps) return false;

  return std::floor(m_t/ft) - m_physFieldFloor > eps;
}

bool
Discretization::fieldrange() const
// *****************************************************************************
//  Decide if physics time falls into a field output time range
//! \return True if physics time falls into a field output time range
// *****************************************************************************
{
  if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;

  const auto eps = std::numeric_limits< tk::real >::epsilon();

  bool output = false;

  const auto& rf = g_inputdeck.get< tag::output, tag::range, tag::field >();
  for (std::size_t i=0; i<rf.size(); ++i)
    if (m_t > rf[i][0] and m_t < rf[i][1])
      output |= std::floor(m_t/rf[i][2]) - m_rangeFieldFloor[i] > eps;

  return output;
}

bool
Discretization::histiter() const
// *****************************************************************************
//  Decide if history output iteration count interval is hit
//! \return True if history output iteration count interval is hit
// *****************************************************************************
{
  const auto hist = g_inputdeck.get< tag::output, tag::iter, tag::history >();
  const auto& hist_points = g_inputdeck.get< tag::history, tag::point >();

  return m_it % hist == 0 and not hist_points.empty();
}

bool
Discretization::histtime() const
// *****************************************************************************
//  Decide if history output physics time interval is hit
//! \return True if history output physics time interval is hit
// *****************************************************************************
{
  if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;

  const auto eps = std::numeric_limits< tk::real >::epsilon();
  const auto ht = g_inputdeck.get< tag::output, tag::time, tag::history >();

  if (ht < eps) return false;

  return std::floor(m_t/ht) - m_physHistFloor > eps;
}

bool
Discretization::histrange() const
// *****************************************************************************
//  Decide if physics time falls into a history output time range
//! \return True if physics time falls into a history output time range
// *****************************************************************************
{
  if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;

  const auto eps = std::numeric_limits< tk::real >::epsilon();

  bool output = false;

  const auto& rh = g_inputdeck.get< tag::output, tag::range, tag::history >();
  for (std::size_t i=0; i<rh.size(); ++i)
    if (m_t > rh[i][0] and m_t < rh[i][1])
      output |= std::floor(m_t/rh[i][2]) - m_rangeHistFloor[i] > eps;

  return output;
}

bool
Discretization::finished() const
// *****************************************************************************
//  Decide if this is the last time step
//! \return True if this is the last time step
// *****************************************************************************
{
  const auto eps = std::numeric_limits< tk::real >::epsilon();
  const auto nstep = g_inputdeck.get< tag::discr, tag::nstep >();
  const auto term = g_inputdeck.get< tag::discr, tag::term >();

  return std::abs(m_t-term) < eps or m_it >= nstep;
}

void
Discretization::status()
// *****************************************************************************
// Output one-liner status report
// *****************************************************************************
{
  // Query after how many time steps user wants TTY dump
  const auto tty = g_inputdeck.get< tag::output, tag::iter, tag::tty >();

  // estimate grind time (taken between this and the previous time step)
  using std::chrono::duration_cast;
  using ms = std::chrono::milliseconds;
  using clock = std::chrono::high_resolution_clock;
  auto grind_time = duration_cast< ms >(clock::now() - m_prevstatus).count();
  m_prevstatus = clock::now();

  if (thisIndex==0 and m_meshid == 0 and not (m_it%tty)) {

    const auto term = g_inputdeck.get< tag::discr, tag::term >();
    const auto t0 = g_inputdeck.get< tag::discr, tag::t0 >();
    const auto nstep = g_inputdeck.get< tag::discr, tag::nstep >();
    const auto diag = g_inputdeck.get< tag::output, tag::iter, tag::diag >();
    const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
    const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
    const auto verbose = g_inputdeck.get< tag::cmd, tag::verbose >();
    const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
    const auto steady = g_inputdeck.get< tag::discr, tag::steady_state >();

    // estimate time elapsed and time for accomplishment
    tk::Timer::Watch ete, eta;
    if (not steady) m_timer.eta( term-t0, m_t-t0, nstep, m_it, ete, eta );

    const auto& def =
      g_inputdeck_defaults.get< tag::cmd, tag::io, tag::screen >();
    tk::Print print( g_inputdeck.get< tag::cmd >().logname( def, m_nrestart ),
                     verbose ? std::cout : std::clog,
                     std::ios_base::app );

    // Output one-liner
    print << std::setfill(' ') << std::setw(8) << m_it << "  "
          << std::scientific << std::setprecision(6)
          << std::setw(12) << m_t << "  "
          << m_dt << "  "
          << std::setfill('0')
          << std::setw(3) << ete.hrs.count() << ":"
          << std::setw(2) << ete.min.count() << ":"
          << std::setw(2) << ete.sec.count() << "  "
          << std::setw(3) << eta.hrs.count() << ":"
          << std::setw(2) << eta.min.count() << ":"
          << std::setw(2) << eta.sec.count() << "  "
          << std::scientific << std::setprecision(6) << std::setfill(' ')
          << std::setw(9) << grind_time << "  ";

    // Augment one-liner status with output indicators
    if (fielditer() or fieldtime() or fieldrange()) print << 'f';
    if (not (m_it % diag)) print << 'd';
    if (histiter() or histtime() or histrange()) print << 't';
    if (m_refined) print << 'h';
    if (not (m_it % lbfreq) && not finished()) print << 'l';
    if (not benchmark && (not (m_it % rsfreq) || finished())) print << 'r';

    if (not m_meshvel_converged) print << 'a';
    m_meshvel_converged = true; // get ready for next time step

    print << std::endl;
  }
}

#include "NoWarning/discretization.def.h"