diff options
| author | Erik Schnetter <schnetter@gmail.com> | 2012-03-18 22:17:29 -0400 |
|---|---|---|
| committer | Barry Wardell <barry.wardell@gmail.com> | 2012-09-11 18:23:12 +0100 |
| commit | c123446d6994268dec4b6b34801ae3267477649f (patch) | |
| tree | e5eeeb3cf30d63d9f4cd1cae1c620b069e41578b /CarpetDev | |
| parent | e03aedc3492ed6ef85468e2795aeb3260033d278 (diff) | |
CarpetIOF5: Continue to implement file reading
Diffstat (limited to 'CarpetDev')
31 files changed, 2363 insertions, 827 deletions
diff --git a/CarpetDev/CarpetIOF5/par/README b/CarpetDev/CarpetIOF5/par/README new file mode 100644 index 000000000..081913573 --- /dev/null +++ b/CarpetDev/CarpetIOF5/par/README @@ -0,0 +1,19 @@ +Example parameter files, listed in order of increasing complexity: + +iof5-uniform.par +iof5-uniform-separate.par +iof5-uniform-fragment.par +iof5-uniform-separate-fragment.par +iof5-refined.par +iof5-refined-cell.par +iof5-multipatch.par +iof5-checkpoint.par +iof5-multipatch-minkowski.par + +iof5-uniform-input.par +iof5-uniform-separate-input.par +iof5-uniform-fragment-input.par +iof5-uniform-separate-fragment-input.par +iof5-refined-input.par +iof5-refined-cell-input.par +iof5-multipatch-input.par diff --git a/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par b/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par index f3bacfa89..583c1dbc6 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par +++ b/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par @@ -47,6 +47,10 @@ CarpetRegrid2::radius_1[2] = 0.2 InitBase::initial_data_setup_method = "init_all_levels" -CarpetIOF5::checkpoint = yes -IO::checkpoint_every = 1024 -IO::checkpoint_dir = $parfile +CarpetIOF5::checkpoint = yes +IO::checkpoint_every = 1024 +IO::checkpoint_ID = yes +IO::checkpoint_on_terminate = yes +IO::checkpoint_dir = $parfile + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch-input.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch-input.par index 202cb179b..d59202db8 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-multipatch-input.par +++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch-input.par @@ -50,3 +50,5 @@ IO::filereader_ID_vars = "grid::coordinates" CarpetIOF5::out_every = 1 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par new file mode 100644 index 000000000..c28ade8a0 --- /dev/null +++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild-input.par @@ -0,0 +1,170 @@ +ActiveThorns = " + ADMBase + ADMCoupling + ADMMacros + Boundary + Carpet + CarpetIOASCII + CarpetIOBasic + CarpetIOF5 + CarpetIOScalar + CarpetLib + CarpetReduce + CarpetRegrid2 + CartGrid3D + CoordBase + CoordGauge + Coordinates + F5 + Formaline + GenericFD + GSL + HDF5 + IDFileADM + InitBase + IOUtil + LoopControl + ML_BSSN + ML_BSSN_Helper + MoL + NewRad + SpaceMask + StaticConformal + SymBase + Time + TimerReport + TmunuBase +" + + + +Cactus::cctk_run_title = "IOF5" +Cactus::cctk_full_warnings = yes +Cactus::highlight_warning_messages = no + +Cactus::cctk_itlast = 0 + +IO::out_dir = $parfile + +Carpet::prolongation_order_space = 5 +Carpet::prolongation_order_time = 2 +driver::ghost_size = 3 +Carpet::use_buffer_zones = yes + +Carpet::domain_from_multipatch = yes +CartGrid3D::type = "multipatch" +Coordinates::coordinate_system = "Thornburg04" +Coordinates::h_cartesian = 0.05 # 0.1 +Coordinates::h_radial = 0.05 # 0.1 +Coordinates::sphere_inner_radius = 1.0 +Coordinates::sphere_outer_radius = 3.0 +Coordinates::n_angular = 40 # 20 +Coordinates::additional_overlap_size = 1 +Coordinates::patch_boundary_size = 3 +Coordinates::outer_boundary_size = 3 + +Carpet::max_refinement_levels = 10 +CarpetRegrid2::num_centres = 1 +CarpetRegrid2::num_levels_1 = 3 +CarpetRegrid2::radius_1[1] = 0.5 +CarpetRegrid2::radius_1[2] = 0.2 + +InitBase::initial_data_setup_method = "init_all_levels" + + + +ADMBase::initial_data = "read from file" +ADMBase::initial_lapse = "read from file" +ADMBase::initial_shift = "read from file" +ADMBase::initial_dtlapse = "read from file" +ADMBase::initial_dtshift = "read from file" + +IO::filereader_ID_dir = "iof5-multipatch-kerrschild" +IO::filereader_ID_files = "iof5-multipatch-kerrschild" +#IO::filereader_ID_vars = " +# ADMBase::metric +# ADMBase::curv +# ADMBase::lapse +# ADMBase::shift +# ADMBase::dtlapse +# ADMBase::dtshift +#" + +ADMBase::evolution_method = "ML_BSSN" +ADMBase::lapse_evolution_method = "ML_BSSN" +ADMBase::shift_evolution_method = "ML_BSSN" +ADMBase::dtlapse_evolution_method = "ML_BSSN" +ADMBase::dtshift_evolution_method = "ML_BSSN" + +ML_BSSN::harmonicN = 1 # 1+log +ML_BSSN::harmonicF = 2.0 # 1+log +ML_BSSN::ShiftGammaCoeff = 0.75 +ML_BSSN::BetaDriver = 1.0 +ML_BSSN::LapseAdvectionCoeff = 1 +ML_BSSN::ShiftAdvectionCoeff = 1 + +ML_BSSN::MinimumLapse = 1.0e-8 + +ML_BSSN::my_initial_boundary_condition = "extrapolate-gammas" +ML_BSSN::my_rhs_boundary_condition = "NewRad" +Boundary::radpower = 2 + +MoL::ODE_Method = "RK4" +MoL::MoL_Intermediate_Steps = 4 +MoL::MoL_Num_Scratch_Levels = 1 + +Time::timestep_method = "given" +Time::timestep = 0.01 + +CarpetIOBasic::outInfo_every = 1 +CarpetIOBasic::outInfo_vars = "ADMBase::lapse" + +CarpetIOScalar::all_reductions_in_one_file = yes +CarpetIOScalar::outScalar_every = 1 +CarpetIOScalar::outScalar_vars = " + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift +" + +CarpetIOASCII::one_file_per_group = yes +CarpetIOASCII::out1D_d = no +CarpetIOASCII::out1D_every = 1024 +CarpetIOASCII::out1D_vars = " + grid::coordinates + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift +" + + + +CarpetIOF5::out_every = 1024 +CarpetIOF5::out_vars = " + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift +" + +CarpetIOF5::checkpoint = yes +IO::checkpoint_ID = yes +IO::checkpoint_every = 1024 +IO::checkpoint_dir = $parfile + + + +Formaline::output_source = no + +TimerReport::out_every = 1024 +TimerReport::out_filename = "TimerReport" +TimerReport::n_top_timers = 50 +TimerReport::output_all_timers_readable = yes diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par new file mode 100644 index 000000000..171ef919e --- /dev/null +++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch-kerrschild.par @@ -0,0 +1,175 @@ +ActiveThorns = " + ADMBase + ADMCoupling + ADMMacros + Boundary + Carpet + CarpetIOASCII + CarpetIOBasic + CarpetIOF5 + CarpetIOScalar + CarpetLib + CarpetReduce + CarpetRegrid2 + CartGrid3D + CoordBase + CoordGauge + Coordinates + Exact + F5 + Formaline + GenericFD + GSL + HDF5 + InitBase + IOUtil + KerrSchild + LoopControl + ML_BSSN + ML_BSSN_Helper + MoL + NewRad + SpaceMask + StaticConformal + SymBase + Time + TimerReport + TmunuBase +" + + + +Cactus::cctk_run_title = "IOF5" +Cactus::cctk_full_warnings = yes +Cactus::highlight_warning_messages = no + +Cactus::cctk_itlast = 0 + +IO::out_dir = $parfile + +Carpet::prolongation_order_space = 5 +Carpet::prolongation_order_time = 2 +driver::ghost_size = 3 +Carpet::use_buffer_zones = yes + +Carpet::domain_from_multipatch = yes +CartGrid3D::type = "multipatch" +Coordinates::coordinate_system = "Thornburg04" +Coordinates::h_cartesian = 0.05 # 0.1 +Coordinates::h_radial = 0.05 # 0.1 +Coordinates::sphere_inner_radius = 1.0 +Coordinates::sphere_outer_radius = 3.0 +Coordinates::n_angular = 40 # 20 +Coordinates::additional_overlap_size = 1 +Coordinates::patch_boundary_size = 3 +Coordinates::outer_boundary_size = 3 + +Carpet::max_refinement_levels = 10 +CarpetRegrid2::num_centres = 1 +CarpetRegrid2::num_levels_1 = 3 +CarpetRegrid2::radius_1[1] = 0.5 +CarpetRegrid2::radius_1[2] = 0.2 + +InitBase::initial_data_setup_method = "init_single_level" +Carpet::init_fill_timelevels = yes + + + +#ADMBase::initial_data = "KerrSchild" +#ADMBase::initial_lapse = "KerrSchild" +#ADMBase::initial_shift = "KerrSchild" +#ADMBase::initial_dtlapse = "KerrSchild" +#ADMBase::initial_dtshift = "KerrSchild" + +#KerrSchild::M = 1.0 +#KerrSchild::a = 0.8 + +ADMBase::initial_data = "Exact" +ADMBase::initial_lapse = "Exact" +ADMBase::initial_shift = "Exact" +ADMBase::initial_dtlapse = "Exact" +ADMBase::initial_dtshift = "Exact" + +Exact::exact_model = "Kerr/Kerr-Schild" +Exact::Kerr_KerrSchild__mass = 1.0 +Exact::Kerr_KerrSchild__spin = 0.8 +Exact::Kerr_KerrSchild__epsilon = 1.0e-4 + +ADMBase::evolution_method = "ML_BSSN" +ADMBase::lapse_evolution_method = "ML_BSSN" +ADMBase::shift_evolution_method = "ML_BSSN" +ADMBase::dtlapse_evolution_method = "ML_BSSN" +ADMBase::dtshift_evolution_method = "ML_BSSN" + +ML_BSSN::harmonicN = 1 # 1+log +ML_BSSN::harmonicF = 2.0 # 1+log +ML_BSSN::ShiftGammaCoeff = 0.75 +ML_BSSN::BetaDriver = 1.0 +ML_BSSN::LapseAdvectionCoeff = 1 +ML_BSSN::ShiftAdvectionCoeff = 1 + +ML_BSSN::MinimumLapse = 1.0e-8 + +ML_BSSN::my_initial_boundary_condition = "extrapolate-gammas" +ML_BSSN::my_rhs_boundary_condition = "NewRad" +Boundary::radpower = 2 + +MoL::ODE_Method = "RK4" +MoL::MoL_Intermediate_Steps = 4 +MoL::MoL_Num_Scratch_Levels = 1 + +Time::timestep_method = "given" +Time::timestep = 0.01 + +CarpetIOBasic::outInfo_every = 1 +CarpetIOBasic::outInfo_vars = "ADMBase::lapse" + +CarpetIOScalar::all_reductions_in_one_file = yes +CarpetIOScalar::outScalar_every = 1 +CarpetIOScalar::outScalar_vars = " + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift +" + +CarpetIOASCII::one_file_per_group = yes +CarpetIOASCII::out1D_d = no +CarpetIOASCII::out1D_every = 1024 +CarpetIOASCII::out1D_vars = " + grid::coordinates + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift +" + + + +CarpetIOF5::out_every = 1024 +CarpetIOF5::out_vars = " + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMBase::dtlapse + ADMBase::dtshift +" + +CarpetIOF5::checkpoint = yes +IO::checkpoint_ID = yes +IO::checkpoint_every = 1024 +IO::checkpoint_dir = $parfile + + + +Formaline::output_source = no + +TimerReport::out_every = 1024 +TimerReport::out_filename = "TimerReport" +TimerReport::n_top_timers = 50 +TimerReport::output_all_timers_readable = yes diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch.par index 9af0afa2f..dbdd1f6b1 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-multipatch.par +++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch.par @@ -46,3 +46,5 @@ InitBase::initial_data_setup_method = "init_all_levels" CarpetIOF5::out_every = 1 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-recover.par b/CarpetDev/CarpetIOF5/par/iof5-recover.par new file mode 100644 index 000000000..ae210ee0f --- /dev/null +++ b/CarpetDev/CarpetIOF5/par/iof5-recover.par @@ -0,0 +1,59 @@ +ActiveThorns = " + Boundary + Carpet + CarpetIOF5 + CarpetLib + CarpetRegrid2 + CartGrid3D + CoordBase + F5 + Formaline + GSL + HDF5 + InitBase + IOUtil + LoopControl + SymBase +" + + + +Cactus::cctk_run_title = "IOF5" +Cactus::cctk_full_warnings = yes +Cactus::highlight_warning_messages = no + +Cactus::cctk_itlast = 2048 + +IO::out_dir = $parfile + +Carpet::domain_from_coordbase = yes +CartGrid3D::type = "coordbase" +CoordBase::domainsize = "minmax" +CoordBase::xmin = -1.0 +CoordBase::ymin = -1.0 +CoordBase::zmin = -1.0 +CoordBase::xmax = +1.0 +CoordBase::ymax = +1.0 +CoordBase::zmax = +1.0 +CoordBase::dx = 0.25 +CoordBase::dy = 0.25 +CoordBase::dz = 0.25 + +Carpet::max_refinement_levels = 10 +CarpetRegrid2::num_centres = 1 +CarpetRegrid2::num_levels_1 = 3 +CarpetRegrid2::radius_1[1] = 0.5 +CarpetRegrid2::radius_1[2] = 0.2 + +InitBase::initial_data_setup_method = "init_all_levels" + +IO::recover = "auto" +IO::recover_dir = "iof5-checkpoint" + +CarpetIOF5::checkpoint = yes +IO::checkpoint_every = 1024 +IO::checkpoint_ID = yes +IO::checkpoint_on_terminate = yes +IO::checkpoint_dir = $parfile + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined-cell-input.par b/CarpetDev/CarpetIOF5/par/iof5-refined-cell-input.par index 322067721..86413e6b1 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-refined-cell-input.par +++ b/CarpetDev/CarpetIOF5/par/iof5-refined-cell-input.par @@ -60,3 +60,5 @@ IO::filereader_ID_vars = "grid::coordinates" CarpetIOF5::out_every = 256 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined-cell.par b/CarpetDev/CarpetIOF5/par/iof5-refined-cell.par index 75dc0c17b..b77536df6 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-refined-cell.par +++ b/CarpetDev/CarpetIOF5/par/iof5-refined-cell.par @@ -56,3 +56,5 @@ InitBase::initial_data_setup_method = "init_all_levels" CarpetIOF5::out_every = 256 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined-input.par b/CarpetDev/CarpetIOF5/par/iof5-refined-input.par index 30917f401..f7e73037d 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-refined-input.par +++ b/CarpetDev/CarpetIOF5/par/iof5-refined-input.par @@ -53,3 +53,5 @@ IO::filereader_ID_vars = "grid::coordinates" CarpetIOF5::out_every = 256 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined.par b/CarpetDev/CarpetIOF5/par/iof5-refined.par index 3fdc219f5..866f0d7f4 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-refined.par +++ b/CarpetDev/CarpetIOF5/par/iof5-refined.par @@ -49,3 +49,5 @@ InitBase::initial_data_setup_method = "init_all_levels" CarpetIOF5::out_every = 256 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment-input.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment-input.par index 3593ea50f..55f546067 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment-input.par +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment-input.par @@ -44,3 +44,5 @@ IO::filereader_ID_vars = "grid::coordinates" CarpetIOF5::out_every = 2 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment.par index 6c363bb52..1e5c9ca19 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment.par +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment.par @@ -42,3 +42,5 @@ CarpetIOF5::fragment_contiguous_components = yes CarpetIOF5::out_every = 2 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-input.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-input.par index 7dfbefd80..3386a062e 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-uniform-input.par +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-input.par @@ -38,9 +38,12 @@ CoordBase::dx = 0.25 CoordBase::dy = 0.25 CoordBase::dz = 0.25 +CarpetIOF5::verbose = yes IO::filereader_ID_dir = "iof5-uniform" IO::filereader_ID_files = "iof5-uniform" IO::filereader_ID_vars = "grid::coordinates" CarpetIOF5::out_every = 2 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment-input.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment-input.par index 55f69a326..e6268ee15 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment-input.par +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment-input.par @@ -44,3 +44,5 @@ IO::filereader_ID_vars = "grid::coordinates" CarpetIOF5::out_every = 2 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment.par index 9f319c9ed..c3efc540c 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment.par +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment.par @@ -43,3 +43,5 @@ CarpetIOF5::fragment_contiguous_components = yes CarpetIOF5::out_every = 2 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-input.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-input.par index 4cfc8ac32..a38ce6f28 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-input.par +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-input.par @@ -44,3 +44,5 @@ IO::filereader_ID_vars = "grid::coordinates" CarpetIOF5::out_every = 2 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate.par index 78d4bd290..39d07475a 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate.par +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate.par @@ -42,3 +42,5 @@ CarpetIOF5::separate_single_component_tensors = yes CarpetIOF5::out_every = 2 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform.par b/CarpetDev/CarpetIOF5/par/iof5-uniform.par index 102096f6b..5779eb157 100644 --- a/CarpetDev/CarpetIOF5/par/iof5-uniform.par +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform.par @@ -40,3 +40,5 @@ CoordBase::dz = 0.25 CarpetIOF5::out_every = 2 CarpetIOF5::out_vars = "grid::coordinates" + +Formaline::output_source = no diff --git a/CarpetDev/CarpetIOF5/param.ccl b/CarpetDev/CarpetIOF5/param.ccl index ee44950de..a002cd016 100644 --- a/CarpetDev/CarpetIOF5/param.ccl +++ b/CarpetDev/CarpetIOF5/param.ccl @@ -4,13 +4,14 @@ SHARES: IO USES STRING filereader_ID_dir AS IO_filereader_ID_dir +USES KEYWORD recover +USES STRING recover_dir AS IO_recover_dir +USES STRING recover_file + USES STRING out_dir AS IO_out_dir USES INT out_every AS IO_out_every - USES INT out_timesteps_per_file -USES STRING recover_dir AS IO_recover_dir - USES BOOLEAN checkpoint_ID USES INT checkpoint_every USES BOOLEAN checkpoint_on_terminate @@ -20,6 +21,12 @@ USES STRING checkpoint_dir AS IO_checkpoint_dir PRIVATE: +BOOLEAN verbose "Verbosity" STEERABLE=always +{ +} "no" + + + STRING filereader_ID_dir "Input directory (overrides IO::filereader_ID_dir)" STEERABLE=always { "^$" :: "Empty: use IO::filereader_ID_dir" @@ -53,17 +60,17 @@ KEYWORD file_content "Create one file for every x" STEERABLE=always "everything" :: "" } "everything" -INT iteration_digits "Minimum number of digits for iteration number" STEERABLE=always -{ - 0:* :: "" -} 10 - STRING out_filename "File name (without extension) for file_content='everything'" STEERABLE=always { "" :: "use the parameter file name" ".+" :: "use this file name" } "" +INT iteration_digits "Minimum number of digits for iteration number" STEERABLE=always +{ + 0:* :: "" +} 8 + INT processor_digits "Minimum number of digits for processor number" STEERABLE=always { 0:* :: "" diff --git a/CarpetDev/CarpetIOF5/schedule.ccl b/CarpetDev/CarpetIOF5/schedule.ccl index 59b86f1cf..da8745b56 100644 --- a/CarpetDev/CarpetIOF5/schedule.ccl +++ b/CarpetDev/CarpetIOF5/schedule.ccl @@ -6,11 +6,19 @@ STORAGE: this_iteration next_output_iteration # Initialisation -SCHEDULE CarpetIOF5_Startup AT startup +SCHEDULE CarpetIOF5_Startup AT startup AFTER Driver_Startup { LANG: C } "Register I/O method" +if (!CCTK_EQUALS(recover, "no")) { + SCHEDULE CarpetIOF5_RecoverParameters AT recover_parameters + { + LANG: C + OPTIONS: level + } "Recover parameters" +} + SCHEDULE CarpetIOF5_Init AT basegrid { LANG: C diff --git a/CarpetDev/CarpetIOF5/src/attributes.cc b/CarpetDev/CarpetIOF5/src/attributes.cc index 9ca1f109c..ff237534a 100644 --- a/CarpetDev/CarpetIOF5/src/attributes.cc +++ b/CarpetDev/CarpetIOF5/src/attributes.cc @@ -20,19 +20,19 @@ namespace CarpetIOF5 { // Write an int attribute - bool WriteAttribute (hid_t const group, - char const *const name, - int const ivalue) + bool WriteAttribute(hid_t const group, + char const *const name, + int const ivalue) { bool error_flag = false; - hid_t const dataspace = FAILWARN (H5Screate (H5S_SCALAR)); + hid_t const dataspace = FAILWARN(H5Screate(H5S_SCALAR)); hid_t const attribute = - FAILWARN (H5Acreate (group, name, H5T_NATIVE_INT, dataspace, - H5P_DEFAULT, H5P_DEFAULT)); - FAILWARN (H5Awrite (attribute, H5T_NATIVE_INT, &ivalue)); - FAILWARN (H5Aclose (attribute)); - FAILWARN (H5Sclose (dataspace)); + FAILWARN(H5Acreate(group, name, H5T_NATIVE_INT, dataspace, + H5P_DEFAULT, H5P_DEFAULT)); + FAILWARN(H5Awrite(attribute, H5T_NATIVE_INT, &ivalue)); + FAILWARN(H5Aclose(attribute)); + FAILWARN(H5Sclose(dataspace)); return error_flag; } @@ -40,19 +40,19 @@ namespace CarpetIOF5 { // Write a double attribute - bool WriteAttribute (hid_t const group, - char const *const name, - double const dvalue) + bool WriteAttribute(hid_t const group, + char const *const name, + double const dvalue) { bool error_flag = false; - hid_t const dataspace = FAILWARN (H5Screate (H5S_SCALAR)); + hid_t const dataspace = FAILWARN(H5Screate(H5S_SCALAR)); hid_t const attribute = - FAILWARN (H5Acreate (group, name, H5T_NATIVE_INT, dataspace, - H5P_DEFAULT, H5P_DEFAULT)); - FAILWARN (H5Awrite (attribute, H5T_NATIVE_DOUBLE, &dvalue)); - FAILWARN (H5Aclose (attribute)); - FAILWARN (H5Sclose (dataspace)); + FAILWARN(H5Acreate(group, name, H5T_NATIVE_INT, dataspace, + H5P_DEFAULT, H5P_DEFAULT)); + FAILWARN(H5Awrite(attribute, H5T_NATIVE_DOUBLE, &dvalue)); + FAILWARN(H5Aclose(attribute)); + FAILWARN(H5Sclose(dataspace)); return error_flag; } @@ -60,44 +60,51 @@ namespace CarpetIOF5 { // Write a string attribute - bool WriteAttribute (hid_t const group, - char const *const name, - char const *const svalue) + bool WriteAttribute(hid_t const group, + char const *const name, + char const *const svalue) { bool error_flag = false; - hid_t const datatype = FAILWARN (H5Tcopy (H5T_C_S1)); - FAILWARN (H5Tset_size (datatype, strlen(svalue) + 1)); - hid_t const dataspace = FAILWARN (H5Screate (H5S_SCALAR)); + hid_t const datatype = FAILWARN(H5Tcopy(H5T_C_S1)); + FAILWARN(H5Tset_size(datatype, strlen(svalue) + 1)); + hid_t const dataspace = FAILWARN(H5Screate(H5S_SCALAR)); hid_t const attribute = - FAILWARN (H5Acreate (group, name, datatype, dataspace, - H5P_DEFAULT, H5P_DEFAULT)); - FAILWARN (H5Awrite (attribute, datatype, svalue)); - FAILWARN (H5Aclose (attribute)); - FAILWARN (H5Sclose (dataspace)); - FAILWARN (H5Tclose (datatype)); + FAILWARN(H5Acreate(group, name, datatype, dataspace, + H5P_DEFAULT, H5P_DEFAULT)); + FAILWARN(H5Awrite(attribute, datatype, svalue)); + FAILWARN(H5Aclose(attribute)); + FAILWARN(H5Sclose(dataspace)); + FAILWARN(H5Tclose(datatype)); return error_flag; } + bool WriteAttribute(hid_t const group, + char const *const name, + string const svalue) + { + return WriteAttribute(group, name, svalue.c_str()); + } + // Write an array of int attributes - bool WriteAttribute (hid_t const group, - char const *const name, - int const *const ivalues, - int const nvalues) + bool WriteAttribute(hid_t const group, + char const *const name, + int const *const ivalues, + int const nvalues) { bool error_flag = false; hsize_t const size = nvalues; - hid_t const dataspace = FAILWARN (H5Screate_simple (1, &size, NULL)); + hid_t const dataspace = FAILWARN(H5Screate_simple(1, &size, NULL)); hid_t const attribute = - FAILWARN (H5Acreate (group, name, H5T_NATIVE_INT, - dataspace, H5P_DEFAULT, H5P_DEFAULT)); - FAILWARN (H5Awrite (attribute, H5T_NATIVE_INT, ivalues)); - FAILWARN (H5Aclose (attribute)); - FAILWARN (H5Sclose (dataspace)); + FAILWARN(H5Acreate(group, name, H5T_NATIVE_INT, + dataspace, H5P_DEFAULT, H5P_DEFAULT)); + FAILWARN(H5Awrite(attribute, H5T_NATIVE_INT, ivalues)); + FAILWARN(H5Aclose(attribute)); + FAILWARN(H5Sclose(dataspace)); return error_flag; } @@ -105,21 +112,21 @@ namespace CarpetIOF5 { // Write an array of double attributes - bool WriteAttribute (hid_t const group, - char const *const name, - double const *const dvalues, - int const nvalues) + bool WriteAttribute(hid_t const group, + char const *const name, + double const *const dvalues, + int const nvalues) { bool error_flag = false; hsize_t const size = nvalues; - hid_t const dataspace = FAILWARN (H5Screate_simple (1, &size, NULL)); + hid_t const dataspace = FAILWARN(H5Screate_simple(1, &size, NULL)); hid_t const attribute = - FAILWARN (H5Acreate (group, name, H5T_NATIVE_DOUBLE, - dataspace, H5P_DEFAULT, H5P_DEFAULT)); - FAILWARN (H5Awrite (attribute, H5T_NATIVE_DOUBLE, dvalues)); - FAILWARN (H5Aclose (attribute)); - FAILWARN (H5Sclose (dataspace)); + FAILWARN(H5Acreate(group, name, H5T_NATIVE_DOUBLE, + dataspace, H5P_DEFAULT, H5P_DEFAULT)); + FAILWARN(H5Awrite(attribute, H5T_NATIVE_DOUBLE, dvalues)); + FAILWARN(H5Aclose(attribute)); + FAILWARN(H5Sclose(dataspace)); return error_flag; } @@ -127,33 +134,33 @@ namespace CarpetIOF5 { // Write an array of string attributes - bool WriteAttribute (hid_t const group, - char const *const name, - char const *const *const svalues, - int const nvalues) + bool WriteAttribute(hid_t const group, + char const *const name, + char const *const *const svalues, + int const nvalues) { bool error_flag = false; size_t maxstrlen = 0; for (int i=0; i<nvalues; ++i) { - maxstrlen = max (maxstrlen, strlen(svalues[i])); + maxstrlen = max(maxstrlen, strlen(svalues[i])); } - vector<char> svalue (nvalues * (maxstrlen+1)); + vector<char> svalue(nvalues * (maxstrlen+1)); for (int i=0; i<nvalues; ++i) { - strncpy (&svalue.at(i*maxstrlen), svalues[i], maxstrlen+1); + strncpy(&svalue.at(i*maxstrlen), svalues[i], maxstrlen+1); } - hid_t const datatype = FAILWARN (H5Tcopy (H5T_C_S1)); - FAILWARN (H5Tset_size (datatype, maxstrlen + 1)); + hid_t const datatype = FAILWARN(H5Tcopy(H5T_C_S1)); + FAILWARN(H5Tset_size(datatype, maxstrlen + 1)); hsize_t const size = nvalues; - hid_t const dataspace = FAILWARN (H5Screate_simple (1, & size, NULL)); + hid_t const dataspace = FAILWARN(H5Screate_simple(1, & size, NULL)); hid_t const attribute = - FAILWARN (H5Acreate (group, name, datatype, dataspace, - H5P_DEFAULT, H5P_DEFAULT)); - FAILWARN (H5Awrite (attribute, datatype, &svalue.front())); - FAILWARN (H5Aclose (attribute)); - FAILWARN (H5Sclose (dataspace)); - FAILWARN (H5Tclose (datatype)); + FAILWARN(H5Acreate(group, name, datatype, dataspace, + H5P_DEFAULT, H5P_DEFAULT)); + FAILWARN(H5Awrite(attribute, datatype, &svalue.front())); + FAILWARN(H5Aclose(attribute)); + FAILWARN(H5Sclose(dataspace)); + FAILWARN(H5Tclose(datatype)); return error_flag; } @@ -161,22 +168,108 @@ namespace CarpetIOF5 { // Write a large string attribute - bool WriteLargeAttribute (hid_t const group, - char const *const name, - char const *const svalue) + bool WriteLargeAttribute(hid_t const group, + char const *const name, + char const *const svalue) { bool error_flag = false; // Create a dataset, since the data may not fit into an attribute - hsize_t const size = strlen (svalue) + 1; - hid_t const dataspace = FAILWARN (H5Screate_simple (1, & size, NULL)); + hsize_t const size = strlen(svalue) + 1; + hid_t const dataspace = FAILWARN(H5Screate_simple(1, & size, NULL)); hid_t const dataset = - FAILWARN (H5Dcreate (group, name, H5T_NATIVE_CHAR, - dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)); - FAILWARN (H5Dwrite (dataset, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL, - H5P_DEFAULT, svalue)); - FAILWARN (H5Dclose (dataset)); - FAILWARN (H5Sclose (dataspace)); + FAILWARN(H5Dcreate(group, name, H5T_NATIVE_CHAR, + dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT)); + FAILWARN(H5Dwrite(dataset, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL, + H5P_DEFAULT, svalue)); + FAILWARN(H5Dclose(dataset)); + FAILWARN(H5Sclose(dataspace)); + + return error_flag; + } + + + + // Read an int attribute + bool ReadAttribute(hid_t const group, + char const * const name, + int& ivalue) + { + bool error_flag = false; + + hid_t const attribute = FAILWARN(H5Aopen_name(group, name)); + FAILWARN(H5Aread(attribute, H5T_NATIVE_INT, &ivalue)); + FAILWARN(H5Aclose(attribute)); + + return error_flag; + } + + + + // Read a double attribute + bool ReadAttribute(hid_t const group, + char const * const name, + double& dvalue) + { + bool error_flag = false; + + hid_t const attribute = FAILWARN(H5Aopen_name(group, name)); + FAILWARN(H5Aread(attribute, H5T_NATIVE_DOUBLE, &dvalue)); + FAILWARN(H5Aclose(attribute)); + + return error_flag; + } + + + + // Read a string attribute + bool ReadAttribute(hid_t const group, + char const * const name, + string& svalue) + { + bool error_flag = false; + + hid_t const attribute = FAILWARN(H5Aopen_name(group, name)); + hid_t const datatype = FAILWARN(H5Aget_type(attribute)); + hsize_t const size = FAILWARN(H5Tget_size(datatype)); + char buf[size+1]; + hid_t const memdatatype = FAILWARN(H5Tcopy(H5T_C_S1)); + FAILWARN(H5Tset_size(memdatatype, size)); + FAILWARN(H5Aread(attribute, memdatatype, buf)); + buf[size] = '\0'; // just to be safe + svalue = string(buf); + FAILWARN(H5Tclose(memdatatype)); + FAILWARN(H5Tclose(datatype)); + FAILWARN(H5Aclose(attribute)); + + return error_flag; + } + + + + // Read a large string attribute + bool ReadLargeAttribute(hid_t const group, + char const *const name, + string& svalue) + { + bool error_flag = false; + + // Read a dataset + hid_t const dataset = FAILWARN(H5Dopen(group, name, H5P_DEFAULT)); + hid_t const dataspace = FAILWARN(H5Dget_space(dataset)); + int const rank = FAILWARN(H5Sget_simple_extent_ndims(dataspace)); + if (rank != 1) error_flag = true; + hsize_t dims[rank]; + FAILWARN(H5Sget_simple_extent_dims(dataspace, dims, NULL)); + hsize_t size = 1; + for (int d=0; d<rank; ++d) size *= dims[d]; + char buf[size+1]; + FAILWARN(H5Dread(dataset, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL, + H5P_DEFAULT, buf)); + buf[size] = '\0'; // just to be safe + svalue = string(buf); + FAILWARN(H5Sclose(dataspace)); + FAILWARN(H5Dclose(dataset)); return error_flag; } diff --git a/CarpetDev/CarpetIOF5/src/distribute.cc b/CarpetDev/CarpetIOF5/src/distribute.cc new file mode 100644 index 000000000..66ad0715f --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/distribute.cc @@ -0,0 +1,627 @@ +#include "distribute.hh" + +#include <cctk_Parameters.h> + +#include <carpet.hh> + +#include <cassert> +#include <cstdlib> +#include <cstring> +#include <vector> + + + +namespace CarpetIOF5 { + + using namespace std; + + + + /*** fragdesc_t *************************************************************/ + + int fragdesc_t::npoints() const + { + int np = 1; + for (int d=0; d<dim; ++d) { + np *= imax[d] - imin[d] + 1; + } + assert(np>0); + return np; + } + + int fragdesc_t::vartypesize() const + { + assert(varindex>=0); + int const vartype = CCTK_VarTypeI(varindex); + assert(vartype>=0); + int const size = CCTK_VarTypeSize(vartype); + assert(size>0); + return size; + } + + MPI_Datatype fragdesc_t::datatype() const + { + assert(varindex>=0); + int const vartype = CCTK_VarTypeI(varindex); + assert(vartype>=0); + switch (vartype) { + case CCTK_VARIABLE_BYTE : return dist::mpi_datatype<CCTK_BYTE >(); + case CCTK_VARIABLE_INT : return dist::mpi_datatype<CCTK_INT >(); + case CCTK_VARIABLE_REAL : return dist::mpi_datatype<CCTK_REAL >(); + case CCTK_VARIABLE_COMPLEX: return dist::mpi_datatype<CCTK_COMPLEX>(); + default: assert(0); + } + } + + + + /*** scatter_t **************************************************************/ + + scatter_t::scatter_t(cGH const *const cctkGH_) + : cctkGH(cctkGH_), + num_received(0), num_sent(0), bytes_allocated(0), + did_send_all(false), + did_receive_sent_all(0), did_receive_all_sent_all(false) + { + DECLARE_CCTK_PARAMETERS; + if (verbose) { + CCTK_INFO("Creating global scatter object"); + } + post_public_recvs(); + } + + + + scatter_t::~scatter_t() + { + DECLARE_CCTK_PARAMETERS; + if (verbose) { + CCTK_INFO("Shutting down global scatter object"); + } + + // Wait until everything has been sent + bool did_send_everything; + for (;;) { + if (verbose) { + CCTK_INFO("Waiting until something has been transmitted..."); + } + did_send_everything = sends.empty(); + if (did_send_everything) break; + do_some_work(true); + } + if (verbose) { + CCTK_INFO("We sent all our data"); + } + + // Notify others: we sent all our data + set_did_send_all(); + + // Wait until everything has been transmitted + bool did_transmit_everything; + for (;;) { + if (verbose) { + CCTK_INFO("Waiting until something has been transmitted..."); + } + did_transmit_everything = + did_receive_all_sent_all and recvs.empty() and sends.empty(); + if (did_transmit_everything) break; + do_some_work(true); + } + if (verbose) { + CCTK_INFO("Everything has been transmitted"); + } + + int const my_difference = num_sent - num_received; + int total_difference; + MPI_Allreduce(const_cast<int*>(&my_difference), &total_difference, + 1, MPI_INT, MPI_SUM, + dist::comm()); + if (total_difference < 0) { + CCTK_WARN(CCTK_WARN_ABORT, + "More received messages than sent messages -- impossible!"); + } + if (total_difference > 0) { + CCTK_WARN(CCTK_WARN_ABORT, "Not all sent messages have been received"); + } + if (verbose) { + CCTK_INFO("Global number of sent and received messages is consistent"); + } + + // Cancel the public receives + if (verbose) { + CCTK_INFO("Cancelling all public receives..."); + } + while (not public_recvs.empty()) { + list<transmission_t*>::iterator const tmi = public_recvs.begin(); + transmission_t *const tm = *tmi; + MPI_Cancel(&tm->request); + public_recvs.erase(tmi); + } + + if (verbose) { + CCTK_INFO("Destroying down global scatter object"); + } + + assert(public_recvs.empty()); + assert(recvs.empty()); + assert(sends.empty()); + } + + + + // Communication tree + int nsiblings() { return 2; } + // Get process id of my first child (the other children are + // consecutive) + int child() { return dist::rank()*nsiblings() + 1; } + // Get process id if my parent + int parent() { return (dist::rank()-1) / nsiblings(); } + + // Check whether we should tell our parent that all our (ours and + // our childrens') messages have been sent + void scatter_t::maybe_send_did_send() + { + DECLARE_CCTK_PARAMETERS; + int need_receive = 0; + for (int p = child(); p<child()+nsiblings(); ++p) { + if (p < dist::size()) ++need_receive; + } + if (did_send_all and did_receive_sent_all == need_receive) { + if (dist::rank() == 0) { + // We are root; now broadcast this to all + set_did_receive_all_sent_all(); + } else { + if (verbose) { + CCTK_INFO("[Telling our parent that we and our children sent all our data]"); + } + to_parent.state = fragdesc_t::state_sent_all; + int const p = parent(); + MPI_Request request; + MPI_Isend(&to_parent, to_parent.num_ints(), MPI_INT, p, tag_desc, + dist::comm(), &request); + } + } + } + + // Broadcast "all message have been sent" to all our children + void scatter_t::send_all_did_send() + { + DECLARE_CCTK_PARAMETERS; + if (verbose) { + CCTK_INFO("[Telling our children that all data have been sent]"); + } + to_children.state = fragdesc_t::state_all_sent_all; + for (int p = child(); p<child()+nsiblings(); ++p) { + if (p < dist::size()) { + MPI_Request request; + MPI_Isend(&to_children, to_children.num_ints(), MPI_INT, p, tag_desc, + dist::comm(), &request); + } + } + } + + // We (this process) sent all our messages + void scatter_t::set_did_send_all() + { + DECLARE_CCTK_PARAMETERS; + if (verbose) { + CCTK_INFO("[We sent all our data]"); + } + assert(not did_send_all); + did_send_all = true; + maybe_send_did_send(); + } + + // One of our children (and all its children) sent all their + // messages + void scatter_t::set_did_receive_sent_all() + { + DECLARE_CCTK_PARAMETERS; + if (verbose) { + CCTK_INFO("[One of our children sent all their data]"); + } + assert(did_receive_sent_all < 2); + ++did_receive_sent_all; + maybe_send_did_send(); + } + + // Our parent broadcast that all messages have been sent + void scatter_t::set_did_receive_all_sent_all() + { + DECLARE_CCTK_PARAMETERS; + if (verbose) { + CCTK_INFO("[Everything has been sent]"); + } + assert(not did_receive_all_sent_all); + did_receive_all_sent_all = true; + send_all_did_send(); + } + + // Dispatch upon a state change message + void scatter_t::handle_state_transition(fragdesc_t const& fd) + { + switch (fd.state) { + case fragdesc_t::state_sent_all: set_did_receive_sent_all(); break; + case fragdesc_t::state_all_sent_all: set_did_receive_all_sent_all(); break; + default: assert(0); + } + } + + + + void scatter_t::send(fragdesc_t const& fd, void const *const data) + { + DECLARE_CCTK_PARAMETERS; + + assert(data); + + if (verbose) { + char *const fullname = CCTK_FullName(fd.varindex); + CCTK_VInfo(CCTK_THORNSTRING, + "Sending data read from variable %s, reflevel %d, map %d, component %d, timelevel %d", + fullname, fd.reflevel, fd.map, fd.src_component, fd.timelevel); + free(fullname); + } + + list<transmission_t*> tosend = split_for_sending(fd, data); + + while (not tosend.empty()) { + list<transmission_t*>::iterator const tmi = tosend.begin(); + transmission_t *const tm = *tmi; + + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, + " Sending to process %d...", tm->fragdesc.process); + } + + // Send descriptor and data + MPI_Request req; + MPI_Isend(&tm->fragdesc, tm->fragdesc.num_ints(), MPI_INT, + tm->fragdesc.process, tag_desc, + dist::comm(), &req); + MPI_Isend(&tm->data[0], tm->fragdesc.npoints(), tm->fragdesc.datatype(), + tm->fragdesc.process, tag_data, + dist::comm(), &tm->request); + + tosend.erase(tmi); + sends.push_back(tm); + } + + // Do some work (if some is available) + do_some_work(); + + if (verbose) { + CCTK_INFO("Done sending"); + } + } + + + + void scatter_t::post_public_recvs() + { + DECLARE_CCTK_PARAMETERS; + + if (verbose) { + CCTK_INFO("Posting public receives"); + } + + // Post public receives + while ((int)public_recvs.size() < num_public_recvs) { + transmission_t *const tm = new transmission_t; + MPI_Irecv(&tm->fragdesc, tm->fragdesc.num_ints(), MPI_INT, + MPI_ANY_SOURCE, tag_desc, + dist::comm(), &tm->request); + public_recvs.push_back(tm); + } + } + + + + void scatter_t::do_some_work(bool const do_wait) + { + DECLARE_CCTK_PARAMETERS; + + if (verbose) { + CCTK_INFO("Checking for progress"); + } + + // Set up an array of all open requests + vector<MPI_Request> requests; + vector<list<transmission_t*>::iterator> iterators; + int const npublic_recvs = public_recvs.size(); + int const nrecvs = recvs.size(); + int const nsends = sends.size(); + int const nrequests = npublic_recvs + nrecvs + nsends; + requests.reserve(nrequests); + iterators.reserve(nrequests); + for (list<transmission_t*>::iterator + tmi = public_recvs.begin(); tmi != public_recvs.end(); ++tmi) + { + transmission_t const *const tm = *tmi; + requests.push_back(tm->request); + iterators.push_back(tmi); + } + for (list<transmission_t*>::iterator + tmi = recvs.begin(); tmi != recvs.end(); ++tmi) + { + transmission_t const *const tm = *tmi; + requests.push_back(tm->request); + iterators.push_back(tmi); + } + for (list<transmission_t*>::iterator + tmi = sends.begin(); tmi != sends.end(); ++tmi) + { + transmission_t const *const tm = *tmi; + requests.push_back(tm->request); + iterators.push_back(tmi); + } + + // Wait for (or test for) some open requests + int outcount; + vector<int> indices(nrequests); + vector<MPI_Status> statuses(nrequests); + if (do_wait) { + if (verbose) { + CCTK_INFO("Waiting for some progress..."); + } + MPI_Waitsome(requests.size(), &requests[0], + &outcount, &indices[0], &statuses[0]); + } else { + if (verbose) { + CCTK_INFO("Testing for some progress..."); + } + MPI_Testsome(requests.size(), &requests[0], + &outcount, &indices[0], &statuses[0]); + } + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, "Completed %d transmissions", outcount); + } + + // Process all completed requests + for (int n=0; n<outcount; ++n) { + int const idx = indices.at(n); + + if (idx < npublic_recvs) { + + // We received a new descriptor + int const source = statuses.at(idx).MPI_SOURCE; + + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, + "Receiving data from process %d", source); + } + + list<transmission_t*>::iterator const tmi = iterators.at(idx); + transmission_t *const tm = *tmi; + public_recvs.erase(tmi); + + if (tm->fragdesc.state != fragdesc_t::state_normal) { + handle_state_transition(tm->fragdesc); + } else { + + // Prepare receiving the data + assert(tm->fragdesc.process == dist::rank()); + tm->data.resize(tm->fragdesc.npoints() * tm->fragdesc.vartypesize()); + bytes_allocated += tm->data.size(); + MPI_Irecv(&tm->data[0], + tm->fragdesc.npoints(), tm->fragdesc.datatype(), + source, tag_data, + dist::comm(), &tm->request); + recvs.push_back(tm); + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, + " Current buffer size: %td bytes", bytes_allocated); + } + + post_public_recvs(); + } + + } else if (idx < npublic_recvs + nrecvs) { + + // We completed receiving a dataset; process it + list<transmission_t*>::iterator const tmi = iterators.at(idx); + transmission_t *const tm = *tmi; + + if (verbose) { + char *const fullname = CCTK_FullName(tm->fragdesc.varindex); + CCTK_VInfo(CCTK_THORNSTRING, + "Completed receiving data for variable %s", fullname); + free(fullname); + } + + write_data(tm); + bytes_allocated -= tm->data.size(); + delete tm; + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, + " Current buffer size: %td bytes", bytes_allocated); + } + recvs.erase(tmi); + ++num_received; + + } else { + + // We completed sending a dataset; forget it + list<transmission_t*>::iterator const tmi = iterators.at(idx); + transmission_t *const tm = *tmi; + + if (verbose) { + char *const fullname = CCTK_FullName(tm->fragdesc.varindex); + CCTK_VInfo(CCTK_THORNSTRING, + "Completed sending data for variable %s", fullname); + free(fullname); + } + + bytes_allocated -= tm->data.size(); + delete tm; + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, + " Current buffer size: %td bytes", bytes_allocated); + } + sends.erase(tmi); + ++num_sent; + + } + } + } + + + + list<scatter_t::transmission_t*> + scatter_t::split_for_sending(fragdesc_t const& fd, void const *const data) + { + DECLARE_CCTK_PARAMETERS; + + int const groupindex = CCTK_GroupIndexFromVarI(fd.varindex); + assert(groupindex>=0); + int const varoffset = fd.varindex - CCTK_FirstVarIndexI(groupindex); + assert(varoffset>=0 and varoffset<=fd.varindex); + + gh const& hh = *Carpet::arrdata.at(groupindex).at(fd.map).hh; + dh const& dd = *Carpet::arrdata.at(groupindex).at(fd.map).dd; + th const& tt = *Carpet::arrdata.at(groupindex).at(fd.map).tt; + + ibbox const& baseext = + hh.baseextents.AT(fd.mglevel).AT(fd.reflevel); + + ibbox const mybox(baseext.lower() + fd.imin * baseext.stride(), + baseext.lower() + fd.imax * baseext.stride(), + baseext.stride()); + + ibset done; + list<transmission_t*> tosend; + dh::light_cboxes const& light_cbox = + dd.light_boxes.at(fd.mglevel).at(fd.reflevel); + for (int c=0; c<hh.components(fd.reflevel); ++c) { + dh::light_dboxes const& light_box = light_cbox.at(c); + ibbox const& intr = light_box.interior; + ibbox const ovlp = mybox & intr; + assert((ovlp & done).empty()); + done += ovlp; + + if (not ovlp.empty()) { + ibbox const& box = ovlp; + transmission_t *const tm = new transmission_t; + tm->fragdesc = fd; + tm->fragdesc.imin = (box.lower() - baseext.lower()) / baseext.stride(); + tm->fragdesc.imax = (box.upper() - baseext.lower()) / baseext.stride(); + tm->fragdesc.component = c; + tm->fragdesc.process = hh.processor(fd.reflevel, c); + + ptrdiff_t const vartypesize = tm->fragdesc.vartypesize(); + ptrdiff_t const ni = tm->fragdesc.imax[0] - tm->fragdesc.imin[0] + 1; + ptrdiff_t const nj = tm->fragdesc.imax[1] - tm->fragdesc.imin[1] + 1; + ptrdiff_t const nk = tm->fragdesc.imax[2] - tm->fragdesc.imin[2] + 1; + ptrdiff_t const di = 1; + ptrdiff_t const dj = ni; + ptrdiff_t const dk = ni * nj; + ptrdiff_t const np = ni * nj * nk; + tm->data.resize(np * vartypesize); + bytes_allocated += tm->data.size(); + if (verbose) { + CCTK_VInfo(CCTK_THORNSTRING, + " Current buffer size: %td bytes", bytes_allocated); + } + char *const dst = &tm->data[0]; + + ptrdiff_t const nis = fd.imax[0] - fd.imin[0] + 1; + ptrdiff_t const njs = fd.imax[1] - fd.imin[1] + 1; + ptrdiff_t const nks = fd.imax[2] - fd.imin[2] + 1; + ptrdiff_t const dis = 1; + ptrdiff_t const djs = nis; + ptrdiff_t const dks = nis * njs; + ptrdiff_t const nps = nis * njs * nks; + ptrdiff_t const i0s = tm->fragdesc.imin[0] - fd.imin[0]; + ptrdiff_t const j0s = tm->fragdesc.imin[1] - fd.imin[1]; + ptrdiff_t const k0s = tm->fragdesc.imin[2] - fd.imin[2]; + ptrdiff_t const ind0s = i0s * dis + j0s * djs + k0s * dks; + char const *const src = &((char const*)data)[vartypesize * ind0s]; + +#pragma omp parallel for collapse(2) + for (ptrdiff_t k=0; k<nk; ++k) { + for (ptrdiff_t j=0; j<nj; ++j) { + ptrdiff_t const ind = j*dj + k*dk; + ptrdiff_t const inds = j*djs + k*dks; + memcpy(&dst[ind], &src[inds], ni*vartypesize); + } + } + + tosend.push_back(tm); + } + } + // Don't enforce this -- mesh refinement boundaries have ghosts + // that do not overlap with any interior + // assert(done == mybox); + + return tosend; + } + + + + void scatter_t::write_data(transmission_t *const tm) + { + fragdesc_t const& fd = tm->fragdesc; + int const groupindex = CCTK_GroupIndexFromVarI(fd.varindex); + assert(groupindex>=0); + int const varoffset = fd.varindex - CCTK_FirstVarIndexI(groupindex); + assert(varoffset>=0 and varoffset<=fd.varindex); + + gh const& hh = *Carpet::arrdata.at(groupindex).at(fd.map).hh; + dh const& dd = *Carpet::arrdata.at(groupindex).at(fd.map).dd; + th const& tt = *Carpet::arrdata.at(groupindex).at(fd.map).tt; + ggf const& ff = + *Carpet::arrdata.at(groupindex).at(fd.map).data.at(varoffset); + int const lc = hh.get_local_component(fd.reflevel, fd.component); + gdata const& data = + *ff.data_pointer(fd.timelevel, fd.reflevel, lc, fd.mglevel); + + ibbox const& baseext = + hh.baseextents.AT(fd.mglevel).AT(fd.reflevel); + + ibbox const mybox(baseext.lower() + fd.imin * baseext.stride(), + baseext.lower() + fd.imax * baseext.stride(), + baseext.stride()); + + dh::light_dboxes const& light_box = + dd.light_boxes.at(fd.mglevel).at(fd.reflevel).at(fd.component); + ibbox const& intr = light_box.interior; + assert(mybox.is_contained_in(intr)); + ibbox const& extr = light_box.exterior; + + ptrdiff_t const vartypesize = fd.vartypesize(); + ptrdiff_t const ni = fd.imax[0] - fd.imin[0] + 1; + ptrdiff_t const nj = fd.imax[1] - fd.imin[1] + 1; + ptrdiff_t const nk = fd.imax[2] - fd.imin[2] + 1; + ptrdiff_t const di = 1; + ptrdiff_t const dj = ni; + ptrdiff_t const dk = ni * nj; + ptrdiff_t const np = ni * nj * nk; + char const *const src = (char const*)&tm->data[0]; + + ivect const lbnd = (extr.lower() - baseext.lower()) / baseext.stride(); + ivect const lsh = extr.shape() / baseext.stride(); + ptrdiff_t const nid = lsh[0]; + ptrdiff_t const njd = lsh[1]; + ptrdiff_t const nkd = lsh[2]; + ptrdiff_t const did = 1; + ptrdiff_t const djd = nid; + ptrdiff_t const dkd = nid * njd; + ptrdiff_t const npd = nid * njd * nkd; + ptrdiff_t const i0d = fd.imin[0] - lbnd[0]; + ptrdiff_t const j0d = fd.imin[1] - lbnd[1]; + ptrdiff_t const k0d = fd.imin[2] - lbnd[2]; + ptrdiff_t const ind0d = i0d * did + j0d * djd + k0d * dkd; + assert(data.has_storage()); + char *const dst = &((char*)data.storage())[ind0d]; + +#pragma omp parallel for collapse(2) + for (ptrdiff_t k=0; k<nk; ++k) { + for (ptrdiff_t j=0; j<nj; ++j) { + ptrdiff_t const indd = j*djd + k*dkd; + ptrdiff_t const ind = j*dj + k*dk; + memcpy(&dst[indd], &src[ind], ni*vartypesize); + } + } + } + +} // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/distribute.hh b/CarpetDev/CarpetIOF5/src/distribute.hh new file mode 100644 index 000000000..420b584c4 --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/distribute.hh @@ -0,0 +1,94 @@ +#ifndef DISTRIBUTE_HH +#define DISTRIBUTE_HH + +#include <cctk.h> + +#include <defs.hh> +#include <vect.hh> + +#include <list> + +#include <mpi.h> + + + +namespace CarpetIOF5 { + + using namespace std; + + + + // Fragment descriptor + struct fragdesc_t { + // only integer entries + int varindex; + static int const mglevel = 0; + int reflevel, map, src_component, timelevel; + ivect imin, imax; // upper bound is inclusive + + // destination + int component, process; + + // meta-messages + enum state_t { state_normal, state_sent_all, state_all_sent_all }; + int state; + + fragdesc_t(): state(state_normal) {} + + int num_ints() const { return sizeof(fragdesc_t) / sizeof(int); } + int npoints() const; + int vartypesize() const; + MPI_Datatype datatype() const; + }; + + + + // Scatter (distribute) Cactus variables + class scatter_t { + + enum tags { tag_desc, tag_data }; + + struct transmission_t { + fragdesc_t fragdesc; + vector<char> data; + MPI_Request request; + }; + + cGH const *const cctkGH; + + // Desired number of public receives to keep posted at all times + static int const num_public_recvs = 10; + list<transmission_t*> public_recvs, recvs, sends; + int num_received, num_sent; + ptrdiff_t bytes_allocated; + + bool did_send_all; + int did_receive_sent_all; + bool did_receive_all_sent_all; + fragdesc_t to_parent, to_children; + + public: + scatter_t(cGH const *const cctkGH_); + ~scatter_t(); + + void send(fragdesc_t const& fd, void const *data); + + private: + void post_public_recvs(); + void do_some_work(bool do_wait = false); + + void maybe_send_did_send(); + void send_all_did_send(); + void set_did_send_all(); + void set_did_receive_sent_all(); + void set_did_receive_all_sent_all(); + void handle_state_transition(fragdesc_t const& fd); + + list<transmission_t*> split_for_sending(fragdesc_t const& fd, + void const *data); + void write_data(transmission_t *const tm); + }; + +} // namespace CarpetIOF5 + +#endif // #ifndef DISTRIBUTE_HH diff --git a/CarpetDev/CarpetIOF5/src/input.cc b/CarpetDev/CarpetIOF5/src/input.cc index da0d7dc00..d42cbd076 100644 --- a/CarpetDev/CarpetIOF5/src/input.cc +++ b/CarpetDev/CarpetIOF5/src/input.cc @@ -10,6 +10,7 @@ #include <limits> #include <sstream> #include <string> +#include <vector> #include <hdf5.h> #include <F5/F5F.h> @@ -30,7 +31,7 @@ // F5 helper namespace { - char const* + char const * ArrayTypeName(ArrayType const array_type) { switch (array_type) { @@ -64,6 +65,8 @@ namespace CarpetIOF5 { cGH const *cctkGH; vector<bool> const input_var; // whether to input this variable + bool const input_past_timelevels; + bool const input_metadata; double time; char const *gridname; @@ -79,11 +82,18 @@ namespace CarpetIOF5 { // string chartname; + scatter_t& scatter; + public: - input_iterator_t (cGH const *const cctkGH_, - vector<bool> const& input_var_) + input_iterator_t(cGH const *const cctkGH_, + vector<bool> const& input_var_, + bool const input_past_timelevels_, + bool const input_metadata_, + scatter_t& scatter_) : cctkGH(cctkGH_), input_var(input_var_), + input_past_timelevels(input_past_timelevels_), + input_metadata(input_metadata_), time(nan), gridname(NULL), topologyname(NULL), index_depth(-1), topological_dimension(-1), @@ -91,7 +101,8 @@ namespace CarpetIOF5 { fieldname(NULL), fieldtype(UnknownArrayType), varindex(-1), fragmentname(NULL), - map(-1), component(-1) + map(-1), component(-1), + scatter(scatter_) { } @@ -99,76 +110,32 @@ namespace CarpetIOF5 { private: - void read_timeslice (F5Path *const path) + void read_timeslice(F5Path *const path) { indent_t indent; cout << indent << "time=" << time << "\n"; -#if 0 - // TODO: Loop over all charts that exist for the grid, or at - // least remember how many maps there are. (This is also written - // as part of the grid structure.) - // F5iterate_grid_atlas? - for (int m=0; m<maps; ++m) { - indent_t indent; - chartname = generate_chartname(cctkGH, m); - cout << indent << "chart=\"" << chartname << "\"\n"; - F5iterate_grids - (path, NULL, grid_iterator, this, NULL, chartname.c_str()); - } -#endif - - F5iterate_grids (path, NULL, grid_iterator, this, NULL, NULL); + F5iterate_grids(path, NULL, grid_iterator, this, NULL, NULL); - // Synchronise - BEGIN_REFLEVEL_LOOP(cctkGH) { -#if 0 - int groups[] = { - CCTK_GroupIndex("GRID::COORDINATES") - }; - int const ngroups = sizeof groups / sizeof *groups; - for (int i=0; i<ngroups; ++i) { - assert (groups[i]>=0); - } - iret = CCTK_SyncGroupsI (cctkGH, ngroups, groups); - assert (iret >= 0); -#endif -#warning "TODO: don't modify boundaries" - BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { - DECLARE_CCTK_ARGUMENTS; -#pragma omp parallel - CCTK_LOOP3_BOUNDARIES(boundaries, cctkGH, - i,j,k, 1,1,1, - cctk_lsh[0]-1,cctk_lsh[1]-1,cctk_lsh[2]-1) - { - int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k); - x[ind3d] = 0.0; - y[ind3d] = 0.0; - z[ind3d] = 0.0; - r[ind3d] = 0.0; - } CCTK_ENDLOOP3_BOUNDARIES(boundaries); - } END_LOCAL_COMPONENT_LOOP; - } END_LOCAL_MAP_LOOP; - } END_REFLEVEL_LOOP; +#warning "TODO: synchronise all read grid functions" } - void read_grid (F5Path *const path) + void read_grid(F5Path *const path) { indent_t indent; - cout << indent << "grid=\"" << gridname << "\"\n"; - // F5iterate_vertex_fields (path, NULL, field_iterator, this, NULL, NULL); - F5iterate_topologies (path, NULL, topology_iterator, this); + cout << indent << "grid=" << gridname << "\n"; + // F5iterate_vertex_fields(path, NULL, field_iterator, this, NULL, NULL); + F5iterate_topologies(path, NULL, topology_iterator, this); } - void read_topology (F5Path *const path) + void read_topology(F5Path *const path) { indent_t indent; herr_t herr; cout << indent - << "topology=\"" << topologyname << "\"" + << "topology=" << topologyname << "" << " (" << (index_depth==0 ? "vertex" : "cell") << ")\n" << indent << "topological dimension=" << topological_dimension << "\n"; @@ -177,13 +144,13 @@ namespace CarpetIOF5 { // Ignore topologies that are only an alias for another topology H5G_stat_t stat; - herr = H5Gget_objinfo (path->Grid_hid, topologyname, false, &stat); - assert (not herr); + herr = H5Gget_objinfo(path->Grid_hid, topologyname, false, &stat); + assert(not herr); if (stat.type == H5G_LINK) { char linkval[100000]; herr = H5Gget_linkval (path->Grid_hid, topologyname, sizeof linkval, linkval); - assert (not herr); + assert(not herr); indent_t indent; cout << indent << "alias for topology \"" << linkval << "\"\n" << indent << "ignoring this topology\n"; @@ -198,16 +165,16 @@ namespace CarpetIOF5 { int const iret = F5LAget_dimensions(path->Topology_hid, FIBER_HDF5_REFINEMENT_INFO, hreffact); - assert (iret == dim); + assert(iret == dim); hsize_t hreffact2[FIBER_MAX_RANK]; void *const pret = F5Tpermute_dimensions(path, dim, hreffact2, hreffact); - assert (pret); + assert(pret); ivect const reffact = h2v(hreffact2); int rl; for (rl=0; rl<reflevels; ++rl) { if (all(reffact == Carpet::spacereffacts.AT(rl))) break; } - assert (rl<reflevels); + assert(rl<reflevels); reflevel = rl; cout << indent << "refinement level " << reflevel << "\n"; @@ -215,13 +182,13 @@ namespace CarpetIOF5 { // F5iterate_topology_fields // (path, NULL, field_iterator, this, chartname.c_str(), NULL); - F5iterate_topology_fields (path, NULL, field_iterator, this, NULL, NULL); + F5iterate_topology_fields(path, NULL, field_iterator, this, NULL, NULL); } - void read_field (F5Path *const path) + void read_field(F5Path *const path) { indent_t indent; - cout << indent << "field=\"" << fieldname << "\"\n"; + cout << indent << "field=" << fieldname << "\n"; interpret_fieldname(cctkGH, fieldname, varindex); #warning "TODO: check all variables in the group" @@ -238,13 +205,14 @@ namespace CarpetIOF5 { assert(iret); // Do we need to iterate over fragments? +#warning "TODO: Should instead check whether attribute FIBER_HDF5_TYPEID_ATTRIB exists (existence indicates fragmentation)" int const is_fragmented = F5Fis_fragmented(path, fieldname); cout << indent << (is_fragmented ? "fragmented" : "not fragmented") << "\n"; if (is_fragmented) { - F5iterate_field_fragments (path, NULL, fragment_iterator, this); + F5iterate_field_fragments(path, NULL, fragment_iterator, this); } else { - read_fragment (path); + read_fragment(path); } F5Fclose(path); @@ -253,17 +221,16 @@ namespace CarpetIOF5 { indent_t indent; cout << indent << "ignoring this field\n"; } +#warning "TODO: keep track of which fields have been read, and complain about unread ones" +#warning "TODO: keep track of which part of a field has been read, and complain about unread parts" } - void read_fragment (F5Path *const path) + void read_fragment(F5Path *const path) { indent_t indent; - int iret; - void *pret; - if (fragmentname) { - cout << indent << "fragment=\"" << fragmentname << "\"\n"; + cout << indent << "fragment=" << fragmentname << "\n"; } else { cout << indent << "no fragment\n"; } @@ -278,7 +245,7 @@ namespace CarpetIOF5 { // Determine map and component from fragment name if (fragmentname) { - interpret_fragmentname (cctkGH, fragmentname, map, component); + interpret_fragmentname(cctkGH, fragmentname, map, component); } else { map = 0; component = CCTK_MyProc(cctkGH); @@ -296,45 +263,51 @@ namespace CarpetIOF5 { hid_t const type_id = F5Fget_type(path); if (F5Tis_convertible(type_id, H5T_NATIVE_DOUBLE)) { cout << indent << "compound type: scalar\n"; - read_variable (path, NULL, varindex); + read_variable(path, NULL, varindex); } else if (F5Tis_convertible(type_id, F5T_VEC3_DOUBLE)) { cout << indent << "compound type: vector\n"; - // This assumes a separated vector; we don't support - // contiguous vectors + // This assumes separated storage; we don't support contiguous + // storage yet assert(is_separated); - read_variable (path, "Dx", varindex+0); - read_variable (path, "Dy", varindex+1); - read_variable (path, "Dz", varindex+2); + read_variable(path, "Dx", varindex+0); + read_variable(path, "Dy", varindex+1); + read_variable(path, "Dz", varindex+2); } else if (F5Tis_convertible(type_id, F5T_METRIC33_DOUBLE)) { cout << indent << "compound type: symmetric tensor\n"; - // Not yet implemented - assert(0); + // This assumes separated storage; we don't support contiguous + // storage yet + assert(is_separated); + read_variable(path, "gxx", varindex+0); + read_variable(path, "gxy", varindex+1); + read_variable(path, "gxz", varindex+2); + read_variable(path, "gyy", varindex+3); + read_variable(path, "gyz", varindex+4); + read_variable(path, "gzz", varindex+5); } else { // Unknown tensor type assert(0); } } - void read_variable (F5Path *const path, char const *const name, - int const var) + void read_variable(F5Path *const path, char const *const name, + int const var) { indent_t indent; herr_t herr; - int ierr; int iret; void *pret; - assert (path); - assert (var >= 0); + assert(path); + assert(var >= 0); - cout << indent << "dataset=\"" << (name ? name : "(null)") << "\"\n"; + cout << indent << "dataset=" << (name ? name : "(null)") << "\n"; - assert (var >= 0); + assert(var >= 0); { char *const fullname = CCTK_FullName(var); - cout << indent << "variable=\"" << fullname << "\"\n"; - free (fullname); + cout << indent << "variable=" << fullname << "\n"; + free(fullname); } @@ -366,16 +339,16 @@ namespace CarpetIOF5 { bool fragment_is_group = false; if (fragmentname) { H5O_info_t info; - herr = H5Oget_info_by_name (path->Field_hid, fragmentname, - &info, H5P_DEFAULT); - assert (not herr); + herr = H5Oget_info_by_name(path->Field_hid, fragmentname, + &info, H5P_DEFAULT); + assert(not herr); fragment_is_group = info.type == H5O_TYPE_GROUP; } cout << indent << "fragment_is_group=" << fragment_is_group << "\n"; hid_t fragment; if (fragment_is_group) { - fragment = H5Gopen (path->Field_hid, fragmentname, H5P_DEFAULT); - assert (fragment >= 0); + fragment = H5Gopen(path->Field_hid, fragmentname, H5P_DEFAULT); + assert(fragment >= 0); } else { fragment = path->Field_hid; } @@ -393,7 +366,7 @@ namespace CarpetIOF5 { H5E_BEGIN_TRY { element = H5Dopen(fragment, datasetname, H5P_DEFAULT); } H5E_END_TRY; - assert (element >= 0); + assert(element >= 0); } bool const field_is_dataset = element < 0; cout << indent << "field_is_dataset=" << field_is_dataset << "\n"; @@ -402,13 +375,13 @@ namespace CarpetIOF5 { // a single element and has already been opened element = fragment; } - assert (element>=0); + assert(element>=0); // Check index depth int index_depth_; iret = F5Tget_index_depth(path, &index_depth_); - assert (iret); - assert (index_depth_ == index_depth); + assert(iret); + assert(index_depth_ == index_depth); // Read the fragment offset. This is stored with the dataset // group. @@ -417,13 +390,13 @@ namespace CarpetIOF5 { hsize_t hoff[FIBER_MAX_RANK]; iret = F5LAget_dimensions(fragment_is_group ? fragment : element, FIBER_FRAGMENT_OFFSET_ATTRIBUTE, hoff); - assert (iret == dim); + assert(iret == dim); hsize_t hoff2[FIBER_MAX_RANK]; pret = F5Tpermute_dimensions(path, dim, hoff2, hoff); - assert (pret); + assert(pret); foff = h2v(hoff2); } - assert (all(foff>=0)); + assert(all(foff>=0)); #if 0 // Read the fragment size. This is stored with the field -- why @@ -432,28 +405,28 @@ namespace CarpetIOF5 { iret = F5LAget_dimensions(path->Field_hid, FIBER_FIELD_DATASPACE_DIMENSIONS_ATTRIBUTE, hlen); - assert (iret == dim); + assert(iret == dim); hsize_t hlen2[FIBER_MAX_RANK]; pret = F5Tpermute_dimensions(path, dim, hlen2, hlen); - assert (pret); + assert(pret); ivect const flen = h2v(hlen2); - assert (all(flen>=0)); + assert(all(flen>=0)); #endif hid_t const space = H5Dget_space(element); - assert (space>=0); + assert(space>=0); iret = H5Sget_simple_extent_ndims(space); - assert (iret == dim); + assert(iret == dim); hsize_t hlen[dim]; iret = H5Sget_simple_extent_dims(space, hlen, NULL); hsize_t hlen2[dim]; pret = F5Tpermute_dimensions(path, dim, hlen2, hlen); - assert (pret); + assert(pret); ivect const flen = h2v(hlen2); - assert (all(flen>=0)); - herr = H5Sclose (space); - assert (not herr); + assert(all(flen>=0)); + herr = H5Sclose(space); + assert(not herr); - ibbox const fbox (foff, foff+flen-1, 1); + ibbox const fbox(foff, foff+flen-1, 1); { indent_t indent; cout << indent @@ -462,131 +435,43 @@ namespace CarpetIOF5 { - // Examine Cactus variable - int const group = CCTK_GroupIndexFromVarI(var); - cGroup groupdata; - ierr = CCTK_GroupData(group, &groupdata); - assert (not ierr); - - // TODO - assert (groupdata.grouptype == CCTK_GF); - assert (groupdata.vartype == CCTK_VARIABLE_REAL); - assert (groupdata.disttype == CCTK_DISTRIB_DEFAULT); - assert (groupdata.stagtype == 0); - assert (groupdata.dim == dim); - - // TODO: This does not work; we need to read the data, and later - // distribute it onto the respective Cactus variables - ENTER_LEVEL_MODE(cctkGH, reflevel) { - ENTER_SINGLEMAP_MODE(cctkGH, map, CCTK_GF) { - ENTER_LOCAL_MODE(cctkGH, component, CCTK_GF) { - DECLARE_CCTK_ARGUMENTS; - - ivect lbnd, lsh, lghosts, ughosts; - ivect imin, imax, ioff, ilen; - for (int d=0; d<dim; ++d) { - lbnd[d] = cctk_lbnd[d]; - lsh[d] = cctk_lsh[d]; - // F5 counts the total overlap, which is the sum of the - // ghost zones on this and on the adjacent component - lghosts[d] = cctk_bbox[2*d ] ? 0 : 2*cctk_nghostzones[d]; - ughosts[d] = cctk_bbox[2*d+1] ? 0 : 2*cctk_nghostzones[d]; - imin[d] = 0; - imax[d] = lsh[d]; - // Do not read in ghost zones - imin[d] += lghosts[d] / 2; - imax[d] -= ughosts[d] / 2; - ioff[d] = lbnd[d] + imin[d]; - ilen[d] = imax[d] - imin[d]; - } - ibbox const lbox (lbnd, lbnd+lsh-1, 1); - ibbox const ibox (ioff, ioff+ilen-1, 1); - ibbox const ovlp = ibox & fbox; - - if (ovlp.empty()) { - indent_t indent; - cout << indent << "dataset does not intersect component " << component << "; skipping component\n"; - continue; - } - - indent_t indent; - cout << indent << "dataset intersects component " << component << "\n" - << indent << "reading bbox " << ovlp.lower() << ":" << ovlp.upper() << " of " << ibox.lower() << ":" << ibox.upper() << "\n"; - - #warning "TODO: assert the whole component is read" - - int const proc = vhh.AT(Carpet::map)->processor(reflevel, component); -#warning "TODO: send dataset to destination process" - assert (proc == dist::rank()); - - int const timelevel = 0; - void *const data = CCTK_VarDataPtrI(cctkGH, timelevel, var); - assert (data); - - - - hvect const hzero(0); - hvect mlsh; - pret = F5Tpermute_dimensions - (path, dim, &mlsh[0], &v2h(lbox.shape())[0]); - assert (pret); - hid_t const memspace = H5Screate_simple(dim, &mlsh[0], NULL); - assert (memspace >= 0); - hvect mmin; - pret = F5Tpermute_dimensions - (path, dim, &mmin[0], &v2h(ovlp.lower()-lbox.lower())[0]); - assert (pret); - hvect mlen; - pret = F5Tpermute_dimensions - (path, dim, &mlen[0], &v2h(ovlp.shape())[0]); - assert (pret); - assert (all(mmin >= hzero)); - assert (all(mmin+mlen <= mlsh)); - herr = H5Sselect_hyperslab - (memspace, H5S_SELECT_SET, &mmin[0], NULL, &mlen[0], NULL); - assert (not herr); - // cout << "mlsh=" << mlsh << " mmin=" << mmin << " mlen=" << mlen << "\n"; - - hvect flsh; - pret = F5Tpermute_dimensions - (path, dim, &flsh[0], &v2h(fbox.shape())[0]); - assert (pret); - hid_t const filespace = H5Screate_simple(dim, &flsh[0], NULL); - assert (filespace >= 0); - hvect fmin; - pret = F5Tpermute_dimensions - (path, dim, &fmin[0], &v2h(ovlp.lower()-fbox.lower())[0]); - assert (pret); - hvect const flen = mlen; - assert (all(fmin >= hzero)); - assert (all(fmin+flen <= flsh)); - herr = H5Sselect_hyperslab - (filespace, H5S_SELECT_SET, &fmin[0], NULL, &flen[0], NULL); - // cout << "flsh=" << flsh << " fmin=" << fmin << " flen=" << flen << "\n"; - - herr = H5Dread (element, H5T_NATIVE_DOUBLE, memspace, filespace, - H5P_DEFAULT, data); - assert (not herr); - - herr = H5Sclose(memspace); - assert (not herr); - herr = H5Sclose(filespace); - assert (not herr); - - } LEAVE_LOCAL_MODE; - } LEAVE_SINGLEMAP_MODE; - } LEAVE_LEVEL_MODE; + fragdesc_t fragdesc; + fragdesc.varindex = var; + fragdesc.reflevel = reflevel; + fragdesc.map = map; + fragdesc.component = component; +#warning "TODO: set timelevel correctly" + fragdesc.timelevel = 0; + fragdesc.imin = fbox.lower(); + fragdesc.imax = fbox.upper(); + + vector<char> data(fragdesc.npoints() * fragdesc.vartypesize()); + + int const vartype = CCTK_VarTypeI(var); + hid_t type; + switch (vartype) { + case CCTK_VARIABLE_INT: type = H5T_NATIVE_INT; break; + case CCTK_VARIABLE_REAL: type = H5T_NATIVE_DOUBLE; break; + default: assert(0); + } + assert(type >= 0); + assert(CCTK_VarTypeSize(vartype) == (int)H5Tget_size(type)); + + herr = H5Dread(element, type, H5S_ALL, H5S_ALL, H5P_DEFAULT, &data[0]); + assert(not herr); + + scatter.send(fragdesc, &data[0]); if (not field_is_dataset) { herr = H5Dclose(element); - assert (not herr); + assert(not herr); } if (fragment_is_group) { herr = H5Gclose(fragment); - assert (not herr); + assert(not herr); } } @@ -594,45 +479,49 @@ namespace CarpetIOF5 { public: - void iterate (hid_t const object) + void iterate(hid_t const object) { - F5iterate_timeslices (object, NULL, timeslice_iterator, this); + F5iterate_timeslices(object, NULL, timeslice_iterator, this); } static - herr_t timeslice_iterator (F5Path *const path, double const time, - void *const userdata) + herr_t timeslice_iterator(F5Path *const path, double const time, + void *const userdata) { input_iterator_t *const iterator = (input_iterator_t*)userdata; iterator->time = time; - iterator->read_timeslice (path); + iterator->read_timeslice(path); + + if (iterator->input_past_timelevels) { + return 0; // continue + } else { + return 1; // done + } return 0; - // finish (successfully) after the first timeslice - // return 1; } static - herr_t grid_iterator (F5Path *const path, char const *const gridname, - void *const userdata) + herr_t grid_iterator(F5Path *const path, char const *const gridname, + void *const userdata) { input_iterator_t *const iterator = (input_iterator_t*)userdata; iterator->gridname = gridname; - iterator->read_grid (path); + iterator->read_grid(path); return 0; } static - herr_t topology_iterator (F5Path *const path, - char const *const topologyname, - int const index_depth, - int const topological_dimension, - void *const userdata) + herr_t topology_iterator(F5Path *const path, + char const *const topologyname, + int const index_depth, + int const topological_dimension, + void *const userdata) { input_iterator_t *const iterator = (input_iterator_t*)userdata; iterator->topologyname = topologyname; iterator->index_depth = index_depth; iterator->topological_dimension = topological_dimension; - iterator->read_topology (path); + iterator->read_topology(path); iterator->topologyname = NULL; iterator->index_depth = -1; iterator->topological_dimension = -1; @@ -640,24 +529,24 @@ namespace CarpetIOF5 { } static - herr_t field_iterator (F5Path *const path, char const *const fieldname, - void *const userdata) + herr_t field_iterator(F5Path *const path, char const *const fieldname, + void *const userdata) { input_iterator_t *const iterator = (input_iterator_t*)userdata; iterator->fieldname = fieldname; - iterator->read_field (path); + iterator->read_field(path); iterator->fieldname = NULL; return 0; } static - herr_t fragment_iterator (F5Path *const path, - char const *const fragmentname, - void *const userdata) + herr_t fragment_iterator(F5Path *const path, + char const *const fragmentname, + void *const userdata) { input_iterator_t *const iterator = (input_iterator_t*)userdata; iterator->fragmentname = fragmentname; - iterator->read_fragment (path); + iterator->read_fragment(path); iterator->fragmentname = NULL; return 0; } @@ -666,12 +555,71 @@ namespace CarpetIOF5 { - void input (cGH const *const cctkGH, - hid_t const file, - vector<bool> const& input_var) + void read_metadata(cGH const *const cctkGH, hid_t const file) + { + DECLARE_CCTK_PARAMETERS; + + assert(cctkGH); + assert(file >= 0); + + herr_t herr; + + + + CCTK_INFO("Reading simulation metadata..."); + + // Open a group holding all metadata + hid_t const group = H5Gopen(file, metadata_group, H5P_DEFAULT); + assert(group >= 0); + + // General information + string fullversion; + ReadAttribute(group, "Cactus version", fullversion); + cout << "Cactus version: " << fullversion << "\n"; + + // Unique identifiers + string config_id; + ReadAttribute(group, "config id", config_id); + cout << "UniqueConfigID: " << config_id << "\n"; + string build_id; + ReadAttribute(group, "build id", build_id); + cout << "UniqueBuildID: " << build_id << "\n"; + string simulation_id; + ReadAttribute(group, "simulation id", simulation_id); + cout << "UniqueSimulationID: " << simulation_id << "\n"; + string run_id; + ReadAttribute(group, "run id", run_id); + cout << "UniqueRunID: " << run_id << "\n"; + + // Parameters + string parameters; + ReadLargeAttribute(group, all_parameters, parameters); + IOUtil_SetAllParameters(parameters.c_str()); + + // Grid structure + string gs; + ReadLargeAttribute(group, grid_structure, gs); +#warning "TODO: set grid structure" + + herr = H5Gclose(group); + assert(not herr); + } + + + + void input(cGH const *const cctkGH, + hid_t const file, + vector<bool> const& input_var, + bool const input_past_timelevels, + bool const input_metadata, + scatter_t& scatter) { - input_iterator_t iterator(cctkGH, input_var); - iterator.iterate (file); +#warning "TODO: not yet implemented" + assert (not input_metadata); + input_iterator_t iterator(cctkGH, input_var, + input_past_timelevels, input_metadata, + scatter); + iterator.iterate(file); } } // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc b/CarpetDev/CarpetIOF5/src/iof5.cc index 388708c6e..01033b24b 100644 --- a/CarpetDev/CarpetIOF5/src/iof5.cc +++ b/CarpetDev/CarpetIOF5/src/iof5.cc @@ -3,6 +3,7 @@ #include <cctk_Parameters.h> #include <util_Table.h> +#include <algorithm> #include <cassert> #include <cstdlib> #include <cstring> @@ -12,12 +13,17 @@ #include <sstream> #include <string> +#ifdef HAVE_DIRENT_H +# include <dirent.h> +#endif + #include <hdf5.h> -#include <iof5.hh> #include "CactusBase/IOUtil/src/ioGH.h" #include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h" +#include "iof5.hh" + namespace CarpetIOF5 { @@ -33,19 +39,19 @@ namespace CarpetIOF5 { // Scheduled startup routine - int CarpetIOF5_Startup () + int CarpetIOF5_Startup() { - CCTK_RegisterBanner ("AMR F5 I/O provided by CarpetIOF5"); + CCTK_RegisterBanner("AMR F5 I/O provided by CarpetIOF5"); - int const GHExtension = CCTK_RegisterGHExtension (CCTK_THORNSTRING); - CCTK_RegisterGHExtensionSetupGH (GHExtension, SetupGH); + int const GHExtension = CCTK_RegisterGHExtension(CCTK_THORNSTRING); + CCTK_RegisterGHExtensionSetupGH(GHExtension, SetupGH); return 0; } // Registered GH extension setup routine - void* SetupGH (tFleshConfig* const fleshconfig, - int const convLevel, cGH* const cctkGH) + void *SetupGH(tFleshConfig *const fleshconfig, + int const convLevel, cGH *const cctkGH) { DECLARE_CCTK_PARAMETERS; @@ -53,18 +59,18 @@ namespace CarpetIOF5 { int const ierr = IOUtil_RegisterRecover("CarpetIOF5 recovery", Input); assert(not ierr); - int const IOMethod = CCTK_RegisterIOMethod ("IOF5"); - CCTK_RegisterIOMethodOutputGH (IOMethod, OutputGH ); - CCTK_RegisterIOMethodTimeToOutput (IOMethod, TimeToOutput ); - CCTK_RegisterIOMethodTriggerOutput (IOMethod, TriggerOutput); - CCTK_RegisterIOMethodOutputVarAs (IOMethod, OutputVarAs ); + int const IOMethod = CCTK_RegisterIOMethod("IOF5"); + CCTK_RegisterIOMethodOutputGH (IOMethod, OutputGH ); + CCTK_RegisterIOMethodTimeToOutput (IOMethod, TimeToOutput ); + CCTK_RegisterIOMethodTriggerOutput(IOMethod, TriggerOutput); + CCTK_RegisterIOMethodOutputVarAs (IOMethod, OutputVarAs ); // there no actual extension data structure return NULL; } // Scheduled initialisation routine - void CarpetIOF5_Init (CCTK_ARGUMENTS) + void CarpetIOF5_Init(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; @@ -76,55 +82,119 @@ namespace CarpetIOF5 { - // Callbacks for CarpetIOHDF5's I/O method - - struct callback_arg_t { - cGH const* cctkGH; + // A mechanism to keep the HDF5 output file open across multiple + // write operations: + hid_t file = H5I_INVALID_HID; + int keep_file_open = 0; + void enter_keep_file_open() + { + ++keep_file_open; + } + void leave_keep_file_open() + { + assert(keep_file_open > 0); + --keep_file_open; + if (keep_file_open==0 and file>=0) { + herr_t const herr = H5Fclose(file); + assert(not herr); + file = H5I_INVALID_HID; + } + } + + + + // Interpret Cactus parameters to decide which variables have been + // selected for output + class selection_t { + int iteration; + int times_set; + vector<bool> selected; + vector<int> last_output_iteration; + static void do_output(int const vindex, + char const *const optstring, + void *const callback_arg) + { + static_cast<selection_t*>(callback_arg)->selected.at(vindex) = true; + } + public: + selection_t(): iteration(-1), times_set(-1) + { + } + bool is_selected(cGH const *const cctkGH, int const vindex) + { + DECLARE_CCTK_PARAMETERS; + // Check whether the parameter out_vars changed only once per + // iteration + if (cctkGH->cctk_iteration != iteration) { + iteration = cctkGH->cctk_iteration; + int const current_times_set = + CCTK_ParameterQueryTimesSet("out_vars", CCTK_THORNSTRING); + // Re-scan out_vars (which is somewhat expensive) only if it + // has changed + if (current_times_set > times_set) { + times_set = current_times_set; + selected.resize(CCTK_NumVars()); + fill(selected.begin(), selected.end(), false); + // for (int n=0; n<CCTK_NumVars(); ++n) selected.at(n) = false; + CCTK_TraverseString(out_vars, do_output, this, CCTK_GROUP_OR_VAR); + } + } + return selected.at(vindex); + } + bool should_output(cGH const *const cctkGH, int const vindex) + { + last_output_iteration.resize(CCTK_NumVars(), -1); + return last_output_iteration.at(vindex) < cctkGH->cctk_iteration; + } + void did_output(cGH const *const cctkGH, int const vindex) + { + last_output_iteration.resize(CCTK_NumVars(), -1); + last_output_iteration.at(vindex) = cctkGH->cctk_iteration; + } }; - void do_output (int const vindex, - char const* const optstring, - void* const callback_arg); + selection_t output_variables; + - int OutputGH (cGH const* const cctkGH) + + // Callbacks for CarpetIOHDF5's I/O method + + int OutputGH(cGH const *const cctkGH) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - static Carpet::Timer timer ("F5::OutputGH"); + static Carpet::Timer timer("F5::OutputGH"); timer.start(); - callback_arg_t callback_arg; - callback_arg.cctkGH = cctkGH; - CCTK_TraverseString (out_vars, do_output, &callback_arg, CCTK_GROUP_OR_VAR); + enter_keep_file_open(); + for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) { + if (TimeToOutput(cctkGH, vindex)) { + TriggerOutput(cctkGH, vindex); + } + } + leave_keep_file_open(); timer.stop(0); return 0; } - void do_output (int const vindex, - char const* const optstring, - void* const callback_arg_) - { - callback_arg_t& callback_arg = - * static_cast<callback_arg_t*>(callback_arg_); - cGH const* const cctkGH = callback_arg.cctkGH; - if (TimeToOutput (cctkGH, vindex)) { - TriggerOutput (cctkGH, vindex); - } - } - - int TimeToOutput (cGH const* const cctkGH, int const vindex) + int TimeToOutput(cGH const *const cctkGH, int const vindex) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; int const numvars = CCTK_NumVars(); - assert (vindex>=0 and vindex<numvars); + assert(vindex>=0 and vindex<numvars); // Output only in global mode if (not do_global_mode) return 0; + // Output only selected variables + if (not output_variables.is_selected(cctkGH, vindex)) return 0; + + if (not output_variables.should_output(cctkGH, vindex)) return 0; + // whether to output at this iteration bool output_this_iteration = false; @@ -144,69 +214,76 @@ namespace CarpetIOF5 { return output_this_iteration ? 1 : 0; } - int TriggerOutput (cGH const* const cctkGH, int const vindex) + int TriggerOutput(cGH const *const cctkGH, int const vindex) { DECLARE_CCTK_PARAMETERS; - char* const fullname = CCTK_FullName(vindex); + char *const fullname = CCTK_FullName(vindex); int const gindex = CCTK_GroupIndexFromVarI(vindex); - char* const groupname = CCTK_GroupName(gindex); - for (char* p=groupname; *p; ++p) *p=tolower(*p); - int const retval = OutputVarAs (cctkGH, fullname, groupname); - free (groupname); - free (fullname); + char *const groupname = CCTK_GroupName(gindex); + for (char *p=groupname; *p; ++p) *p=tolower(*p); + int const retval = OutputVarAs(cctkGH, fullname, groupname); + free(groupname); + free(fullname); + + output_variables.did_output(cctkGH, vindex); return retval; } - int OutputVarAs (cGH const* const cctkGH, - const char* const varname, const char* const alias) + int OutputVarAs(cGH const *const cctkGH, + const char *const varname, const char *const alias) { DECLARE_CCTK_PARAMETERS; - assert (is_level_mode()); + assert(is_level_mode()); BEGIN_GLOBAL_MODE(cctkGH) { DECLARE_CCTK_ARGUMENTS; - CCTK_VInfo (CCTK_THORNSTRING, - "F5::OutputVarAs: iteration=%d", cctk_iteration); + CCTK_VInfo(CCTK_THORNSTRING, + "F5::OutputVarAs: iteration=%d, variable=%s", + cctk_iteration, varname); // We don't know how to open multiple files yet - assert (CCTK_EQUALS (file_content, "everything")); + assert(CCTK_EQUALS(file_content, "everything")); // Open file static bool first_time = true; // The file name doesn't matter since we currently write // everything into a single file - int const vindex = CCTK_VarIndex (varname); - assert (vindex >= 0); - string const basename = generate_basename (cctkGH, vindex); + int const vindex = CCTK_VarIndex(varname); + assert(vindex >= 0); + string const basename = generate_basename(cctkGH, vindex); int const myproc = CCTK_MyProc(cctkGH); int const proc = myproc; string const name = - create_filename (cctkGH, basename, proc, io_dir_output, first_time); + create_filename(cctkGH, basename, cctkGH->cctk_iteration, proc, + io_dir_output, first_time); indent_t indent; cout << indent << "process=" << proc << "\n"; + enter_keep_file_open(); bool const truncate_file = first_time and IO_TruncateOutputFiles(cctkGH); - hid_t const file = - truncate_file ? - H5Fcreate (name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) : - H5Fopen (name.c_str(), H5F_ACC_RDWR , H5P_DEFAULT); - assert (file >= 0); + if (file < 0) { + // Reuse file hid if file is already open + file = + truncate_file ? + H5Fcreate(name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) : + H5Fopen (name.c_str(), H5F_ACC_RDWR , H5P_DEFAULT); + assert(file >= 0); + } first_time = false; vector<bool> output_var(CCTK_NumVars()); output_var.at(vindex) = true; - output (cctkGH, file, output_var, false); + output(cctkGH, file, output_var, false, true); // Close file - herr_t const herr = H5Fclose (file); - assert (not herr); + leave_keep_file_open(); } END_GLOBAL_MODE; @@ -217,83 +294,95 @@ namespace CarpetIOF5 { // Checkpointing - void Checkpoint (cGH const* const cctkGH, int const called_from) + void Checkpoint(cGH const *const cctkGH, int const called_from) { - assert (is_global_mode()); + DECLARE_CCTK_PARAMETERS; + + assert(is_global_mode()); + + CCTK_VInfo(CCTK_THORNSTRING, + "F5::Checkpoint: iteration=%d", cctkGH->cctk_iteration); #if 0 // generate filenames for both the temporary and real checkpoint // files int const ioproc = CCTK_MyProc(cctkGH); int const parallel_io = 1; - char* const filename = - IOUtil_AssembleFilename (cctkGH, NULL, "", ".f5", - called_from, ioproc, not parallel_io); - char* const tempname = - IOUtil_AssembleFilename (cctkGH, NULL, ".tmp", ".f5", - called_from, ioproc, not parallel_io); + char *const filename = + IOUtil_AssembleFilename(cctkGH, NULL, "", ".f5", + called_from, ioproc, not parallel_io); + char *const tempname = + IOUtil_AssembleFilename(cctkGH, NULL, ".tmp", ".f5", + called_from, ioproc, not parallel_io); #endif int const myproc = CCTK_MyProc(cctkGH); int const proc = myproc; string const name = - create_filename (cctkGH, "checkpoint", proc, io_dir_checkpoint, true); + create_filename(cctkGH, "checkpoint", cctkGH->cctk_iteration, proc, + io_dir_checkpoint, true); string const tempname = - create_filename (cctkGH, "checkpoint.tmp", proc, io_dir_checkpoint, true); + create_filename(cctkGH, "checkpoint.tmp", cctkGH->cctk_iteration, proc, + io_dir_checkpoint, true); hid_t const file = - H5Fcreate (tempname.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - assert (file >= 0); + H5Fcreate(tempname.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + assert(file >= 0); - vector<bool> output_var(CCTK_NumVars()); + vector<bool> output_var(CCTK_NumVars(), false); for (int gindex=0; gindex<CCTK_NumGroups(); ++gindex) { - // only checkpoint groups with storage - if (CCTK_QueryGroupStorageI(cctkGH, gindex) <= 0) continue; - if (CCTK_NumVarsInGroupI(gindex) == 0) continue; + // We can't check for storage here, since this requires at least + // level mode + if (not CCTK_QueryGroupStorageI(cctkGH, gindex)) continue; // do not checkpoint groups with a "checkpoint=no" tag cGroup gdata; - CCTK_GroupData (gindex, &gdata); + CCTK_GroupData(gindex, &gdata); int const len = - Util_TableGetString (gdata.tagstable, 0, NULL, "checkpoint"); + Util_TableGetString(gdata.tagstable, 0, NULL, "checkpoint"); if (len > 0) { char buf[1000]; - Util_TableGetString (gdata.tagstable, sizeof buf, buf, "checkpoint"); + Util_TableGetString(gdata.tagstable, sizeof buf, buf, "checkpoint"); if (CCTK_EQUALS(buf, "no")) continue; - assert (CCTK_EQUALS(buf, "yes")); + assert(CCTK_EQUALS(buf, "yes")); } - int const first_vindex = CCTK_FirstVarIndexI (gindex); - int const num_vars = CCTK_NumVarsInGroupI (gindex); - for (int vindex=first_vindex; vindex<first_vindex+num_vars; ++vindex) { - output_var.at(vindex) = true; + int const first_vindex = CCTK_FirstVarIndexI(gindex); + int const num_vars = CCTK_NumVarsInGroupI(gindex); + if (num_vars > 0) { + for (int vindex=first_vindex; vindex<first_vindex+num_vars; ++vindex) { + output_var.at(vindex) = true; + } } } - output (cctkGH, file, output_var, true); + // NOTE: We could write the metadata into only one of the process + // files (to save space), or write it into a separate metadata + // file + output(cctkGH, file, output_var, true, true); // Close file - herr_t const herr = H5Fclose (file); - assert (not herr); + herr_t const herr = H5Fclose(file); + assert(not herr); // Wait until all files have been written, then rename the // checkpoint files // TODO: ensure there were no errors - CCTK_Barrier (cctkGH); - int const ierr = rename (tempname.c_str(), name.c_str()); - assert (not ierr); + CCTK_Barrier(cctkGH); + int const ierr = rename(tempname.c_str(), name.c_str()); + assert(not ierr); } - void CarpetIOF5_InitialDataCheckpoint (CCTK_ARGUMENTS) + void CarpetIOF5_InitialDataCheckpoint(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; if (checkpoint and checkpoint_ID) { - Checkpoint (cctkGH, CP_INITIAL_DATA); + Checkpoint(cctkGH, CP_INITIAL_DATA); } } - void CarpetIOF5_EvolutionCheckpoint (CCTK_ARGUMENTS) + void CarpetIOF5_EvolutionCheckpoint(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; @@ -305,11 +394,11 @@ namespace CarpetIOF5 { checkpoint and (checkpoint_by_iteration or checkpoint_next); if (do_checkpoint) { - Checkpoint (cctkGH, CP_EVOLUTION_DATA); + Checkpoint(cctkGH, CP_EVOLUTION_DATA); } } - void CarpetIOF5_TerminationCheckpoint (CCTK_ARGUMENTS) + void CarpetIOF5_TerminationCheckpoint(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; @@ -320,15 +409,18 @@ namespace CarpetIOF5 { bool const do_checkpoint = not did_checkpoint; if (do_checkpoint) { - Checkpoint (cctkGH, CP_EVOLUTION_DATA); + Checkpoint(cctkGH, CP_EVOLUTION_DATA); } } } - int Input (cGH* const cctkGH, - char const* const basefilename, int const called_from) + // Recovery information + static int recovery_iteration = -1; + + int Input(cGH *const cctkGH, + char const *const basefilename, int const called_from) { DECLARE_CCTK_PARAMETERS; @@ -336,27 +428,31 @@ namespace CarpetIOF5 { - assert (is_level_mode()); BEGIN_GLOBAL_MODE(cctkGH) { DECLARE_CCTK_ARGUMENTS; - CCTK_VInfo (CCTK_THORNSTRING, "F5::Input: iteration=%d", cctk_iteration); - - - - assert (called_from == CP_RECOVER_PARAMETERS or - called_from == CP_RECOVER_DATA or - called_from == FILEREADER_DATA); + assert(called_from == CP_RECOVER_PARAMETERS or + called_from == CP_RECOVER_DATA or + called_from == FILEREADER_DATA); bool const in_recovery = called_from == CP_RECOVER_PARAMETERS or called_from == CP_RECOVER_DATA; - // We don't know how to do this yet - assert (called_from != CP_RECOVER_PARAMETERS); + if (in_recovery) { + CCTK_VInfo(CCTK_THORNSTRING, "F5::Input: recovering iteration %d", + recovery_iteration); + } else { + CCTK_VInfo(CCTK_THORNSTRING, "F5::Input: reading iteration %d", + cctk_iteration); + } + cout << "called_from=" << + (called_from == CP_RECOVER_PARAMETERS ? "CP_RECOVER_PARAMETERS" : + called_from == CP_RECOVER_DATA ? "CP_RECOVER_DATA" : + called_from == FILEREADER_DATA ? "FILEREADER_DATA" : + NULL) << "\n"; // Determine which variables to read - ioGH const* const ioUtilGH = - (ioGH const*) CCTK_GHExtension (cctkGH, "IO"); + ioGH const *const ioUtilGH = (ioGH const*)CCTK_GHExtension(cctkGH, "IO"); vector<bool> input_var(CCTK_NumVars(), true); if (ioUtilGH->do_inVars) { for (int n=0; n<CCTK_NumVars(); ++n) { @@ -366,11 +462,13 @@ namespace CarpetIOF5 { + scatter_t scatter(cctkGH); + // Open file // string const basename = // in_recovery // ? "checkpoint" - // : generate_basename (cctkGH, CCTK_VarIndex("grid::r")); + // : generate_basename(cctkGH, CCTK_VarIndex("grid::r")); string const basename = basefilename; // Keep track of which files could be read, and which could not @@ -382,8 +480,8 @@ namespace CarpetIOF5 { // Loop over all (possible) files for (int proc=myproc; ; proc+=nprocs) { string const name = - create_filename (cctkGH, basename, proc, - in_recovery ? io_dir_recover : io_dir_input, false); + create_filename(cctkGH, basename, cctkGH->cctk_iteration, proc, + in_recovery ? io_dir_recover : io_dir_input, false); bool file_exists; H5E_BEGIN_TRY { @@ -398,15 +496,19 @@ namespace CarpetIOF5 { indent_t indent; cout << indent << "process=" << proc << "\n"; - hid_t const file = H5Fopen (name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); - assert (file >= 0); + hid_t const file = H5Fopen(name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + assert(file >= 0); // Iterate over all time slices - input (cctkGH, file, input_var); + bool const input_past_timelevels = in_recovery; +#warning "TODO: read metadata when recoverying parameters" + bool const input_metadata = false; + input(cctkGH, file, input_var, input_past_timelevels, input_metadata, + scatter); // Close file - herr = H5Fclose (file); - assert (not herr); + herr = H5Fclose(file); + assert(not herr); } { @@ -415,9 +517,10 @@ namespace CarpetIOF5 { dist::comm()); if (maxfoundproc == -1) { string const name = - create_filename (cctkGH, basename, notfoundproc, - in_recovery ? io_dir_recover : io_dir_input, - false); + create_filename(cctkGH, basename, + cctkGH->cctk_iteration, notfoundproc, + in_recovery ? io_dir_recover : io_dir_input, + false); CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, "Could not read input file \"%s\"", name.c_str()); return 1; @@ -434,4 +537,91 @@ namespace CarpetIOF5 { return 0; // no error } + int CarpetIOF5_RecoverParameters() + { +#if 0 + return IOUtil_RecoverParameters(Input, ".f5", "F5"); +#endif + + DECLARE_CCTK_PARAMETERS; + + char const *const IO_dir = IO_recover_dir; + char const *const F5_dir = recover_dir; + bool const use_IO_dir = strcmp(F5_dir, "") == 0; + char const *const my_recover_dir = use_IO_dir ? IO_dir : F5_dir; + + DIR *const dir = opendir(my_recover_dir); + if (not dir) { + // The recovery directory does not exist + if (CCTK_Equals(recover, "autoprobe")) { + // This is harmless when "autoprobe" is used + CCTK_VInfo(CCTK_THORNSTRING, + "Recovery directory \"%s\" doesn't exist", my_recover_dir); + return 0; + } else { + // This is an error when "auto" is used + CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Recovery directory \"%s\" doesn't exist", my_recover_dir); + return -2; + } + } + + // Get the list of potential recovery files + char const *const my_recover_file = "checkpoint"; + string const prefix = string(my_recover_file) + ".i"; + string const infix = ".p"; + string const suffix = ".f5"; + assert(recovery_iteration < 0); + while (dirent *const file = readdir(dir)) { + char *p = file->d_name; + + // First check the file prefix + if (prefix.compare(0, prefix.length(), p, prefix.length()) != 0) continue; + p += prefix.length(); + + // Now check if there is an iteration number following the file + // prefix + int const iter = strtol(p, &p, 10); + if (!*p) continue; + + // Read the process number + if (infix.compare(0, infix.length(), p, infix.length()) != 0) continue; + p += infix.length(); + int const proc = strtol(p, &p, 10); + if (!*p) continue; + + // Finally check the file extension + if (suffix.compare(0, suffix.length(), p, suffix.length()) != 0) continue; + p += suffix.length(); + + // Check whether we read the whole string + if (*p) continue; + + // Found a recovery file by that basename + recovery_iteration = max(recovery_iteration, iter); + } + closedir(dir); + + // There is no recovery file + if (recovery_iteration < 0) { + if (CCTK_Equals(recover, "autoprobe")) { + // This is harmless when "autoprobe" is used + CCTK_VInfo(CCTK_THORNSTRING, + "No F5 checkpoint files with basefilename \"%s\" found in " + "recovery directory \"%s\"", + my_recover_file, my_recover_dir); + return 0; + } else { + // This is an error when "auto" is used + CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING, + "No F5 checkpoint files with basefilename \"%s\" found in " + "recovery directory \"%s\"", + my_recover_file, my_recover_dir); + return -1; + } + } + + return Input(NULL, "checkpoint", CP_RECOVER_PARAMETERS); + } + } // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/iof5.hh b/CarpetDev/CarpetIOF5/src/iof5.hh index 7e54eb9c9..352c83224 100644 --- a/CarpetDev/CarpetIOF5/src/iof5.hh +++ b/CarpetDev/CarpetIOF5/src/iof5.hh @@ -24,6 +24,8 @@ #include <carpet.hh> +#include "distribute.hh" + // This requires defining a boolean local variable "error_flag", @@ -33,13 +35,13 @@ template<typename T> static -T failwarn (bool& error_flag, T const expr, - int const line, char const* const file, char const* const thorn, - char const* const msg) +T failwarn(bool& error_flag, T const expr, + int const line, char const *const file, char const *const thorn, + char const *const msg) { if (expr < 0) { - CCTK_VWarn (CCTK_WARN_ALERT, line, file, thorn, - "Expression \"%s\" return %d", msg, (int)expr); + CCTK_VWarn(CCTK_WARN_ALERT, line, file, thorn, + "Expression \"%s\" return %d", msg, (int)expr); error_flag = true; } return expr; @@ -84,9 +86,9 @@ namespace CarpetIOF5 { CCTK_REAL const nan = numeric_limits<CCTK_REAL>::quiet_NaN(); // Special group and attribute names - char const* const metadata_group = "Parameters and Global Attributes"; - char const* const all_parameters = "All Parameters"; - extern char const* const grid_structure; + char const *const metadata_group = "Parameters and Global Attributes"; + char const *const all_parameters = "All Parameters"; + extern char const *const grid_structure; @@ -102,17 +104,17 @@ namespace CarpetIOF5 { // Conversion operators for these datatypes static inline - hvect v2h (ivect const& a) + hvect v2h(ivect const& a) { return hvect(a); } static inline - ivect h2v (hvect const& a) + ivect h2v(hvect const& a) { return ivect(a); } static inline - ivect h2v (hsize_t const* const a) + ivect h2v(hsize_t const *const a) { ivect res; for (int d=0; d<dim; ++d) { @@ -122,7 +124,7 @@ namespace CarpetIOF5 { } static inline - F5_vec3_point_t v2p (rvect const& a) + F5_vec3_point_t v2p(rvect const& a) { F5_vec3_point_t res; res.x = a[0]; @@ -132,7 +134,7 @@ namespace CarpetIOF5 { } static inline - F5_vec3_float_t v2f (rvect const& a) + F5_vec3_float_t v2f(rvect const& a) { F5_vec3_float_t res; res.x = a[0]; @@ -142,7 +144,7 @@ namespace CarpetIOF5 { } static inline - F5_vec3_double_t v2d (rvect const& a) + F5_vec3_double_t v2d(rvect const& a) { F5_vec3_double_t res; res.x = a[0]; @@ -154,120 +156,142 @@ namespace CarpetIOF5 { // Handle HDF5 attributes more comfortably - bool WriteAttribute (hid_t const group, - char const* const name, - int const ivalue); - bool WriteAttribute (hid_t const group, - char const* const name, - double const dvalue); - bool WriteAttribute (hid_t const group, - char const* const name, - char const* const svalue); - bool WriteAttribute (hid_t const group, - char const* const name, - int const* const ivalues, - int const nvalues); - bool WriteAttribute (hid_t const group, - char const* const name, - double const* const dvalues, - int const nvalues); - bool WriteAttribute (hid_t const group, - char const* const name, - char const* const* const svalues, - int const nvalues); - bool WriteLargeAttribute (hid_t const group, - char const* const name, - char const* const svalue); + bool WriteAttribute(hid_t const group, + char const * const name, + int const ivalue); + bool WriteAttribute(hid_t const group, + char const * const name, + double const dvalue); + bool WriteAttribute(hid_t const group, + char const *const name, + char const *const svalue); + bool WriteAttribute(hid_t const group, + char const *const name, + string const svalue); + bool WriteAttribute(hid_t const group, + char const *const name, + int const *const ivalues, + int const nvalues); + bool WriteAttribute(hid_t const group, + char const *const name, + double const *const dvalues, + int const nvalues); + bool WriteAttribute(hid_t const group, + char const *const name, + char const *const *const svalues, + int const nvalues); + bool WriteLargeAttribute(hid_t const group, + char const *const name, + char const *const svalue); + + bool ReadAttribute(hid_t const group, + char const * const name, + int& ivalue); + bool ReadAttribute(hid_t const group, + char const * const name, + double& dvalue); + bool ReadAttribute(hid_t const group, + char const *const name, + string& svalue); + bool ReadLargeAttribute(hid_t const group, + char const *const name, + string& svalue); // Generate a good file name ("alias") for a variable string - generate_basename (cGH const* const cctkGH, - int const variable); + generate_basename(cGH const *const cctkGH, + int const variable); // Create the final file name on a particular processor enum io_dir_t {io_dir_input, io_dir_output, io_dir_recover, io_dir_checkpoint}; string - create_filename (cGH const* const cctkGH, - string const basename, - int const proc, - io_dir_t const io_dir, - bool const create_directories); + create_filename(cGH const *const cctkGH, + string const basename, + int const iteration, + int const proc, + io_dir_t const io_dir, + bool const create_directories); // Generate a good grid name (simulation name) string - generate_gridname (cGH const* const cctkGH); + generate_gridname(cGH const *const cctkGH); // Generate a good topology name (refinement level name) string - generate_topologyname (cGH const* const cctkGH, - int const gi, - ivect const& reffact); + generate_topologyname(cGH const *const cctkGH, + int const gi, + ivect const& reffact); // Generate a good chart name (tensor basis name) string - generate_chartname (cGH const* const cctkGH); + generate_chartname(cGH const *const cctkGH); // Generate a good fragment name (map and component name) string - generate_fragmentname (cGH const* const cctkGH, int const m, int const c); + generate_fragmentname(cGH const *const cctkGH, int const m, int const c); void - interpret_fragmentname (cGH const* const cctkGH, - char const* const fragmentname, - int& m, int& c); + interpret_fragmentname(cGH const *const cctkGH, + char const *const fragmentname, + int& m, int& c); // Generate a good field name (group or variable name) string - generate_fieldname (cGH const* const cctkGH, - int const vi, tensortype_t const tt); + generate_fieldname(cGH const *const cctkGH, + int const vi, tensortype_t const tt); void interpret_fieldname(cGH const *const cctkGH, string fieldname, int& vi); void - write_metadata (cGH const* const cctkGH, hid_t const file); + write_metadata(cGH const *const cctkGH, hid_t const file); // Handle Carpet's grid structure (this should move to Carpet and/or // CarpetLib) string - serialise_grid_structure (cGH const* const cctkGH); - + serialise_grid_structure(cGH const *const cctkGH); - void output (cGH const* const cctkGH, - hid_t const file, - vector<bool> const& output_var, - bool const output_everything); - void input (cGH const* const cctkGH, + void output(cGH const *const cctkGH, hid_t const file, - vector<bool> const& input_var); + vector<bool> const& output_var, + bool const output_past_timelevels, + bool const output_metadata); + + void input(cGH const *const cctkGH, + hid_t const file, + vector<bool> const& input_var, + bool const input_past_timelevels, + bool const input_metadata, + scatter_t& scatter); // Scheduled routines extern "C" { - int CarpetIOF5_Startup (); - void CarpetIOF5_Init (CCTK_ARGUMENTS); - void CarpetIOF5_InitialDataCheckpoint (CCTK_ARGUMENTS); - void CarpetIOF5_EvolutionCheckpoint (CCTK_ARGUMENTS); - void CarpetIOF5_TerminationCheckpoint (CCTK_ARGUMENTS); + int CarpetIOF5_Startup(); + int CarpetIOF5_RecoverParameters(); + void CarpetIOF5_Init(CCTK_ARGUMENTS); + void CarpetIOF5_InitialDataCheckpoint(CCTK_ARGUMENTS); + void CarpetIOF5_EvolutionCheckpoint(CCTK_ARGUMENTS); + void CarpetIOF5_TerminationCheckpoint(CCTK_ARGUMENTS); } // Registered GH extension setup routine - void* SetupGH (tFleshConfig* const fleshconfig, - int const convLevel, cGH* const cctkGH); + void *SetupGH(tFleshConfig *const fleshconfig, + int const convLevel, cGH *const cctkGH); // Callbacks for CarpetIOF5's I/O method - int Input (cGH* const cctkGH, - char const* const basefilename, int const called_from); - int OutputGH (cGH const* const cctkGH); - int TimeToOutput (cGH const* const cctkGH, int const vindex); - int TriggerOutput (cGH const* const cctkGH, int const vindex); - int OutputVarAs (cGH const* const cctkGH, - const char* const varname, const char* const alias); + int Input(cGH *const cctkGH, + char const *const basefilename, int const called_from); + int OutputGH(cGH const *const cctkGH); + int TimeToOutput(cGH const *const cctkGH, int const vindex); + int TriggerOutput(cGH const *const cctkGH, int const vindex); + int OutputVarAs(cGH const *const cctkGH, + const char *const varname, const char *const alias); } // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/make.code.defn b/CarpetDev/CarpetIOF5/src/make.code.defn index 4603ab0bc..8bb25d479 100644 --- a/CarpetDev/CarpetIOF5/src/make.code.defn +++ b/CarpetDev/CarpetIOF5/src/make.code.defn @@ -1,7 +1,7 @@ # Main make.code.defn file for thorn CarpetIOF5 # Source files in this directory -SRCS = attributes.cc input.cc iof5.cc output.cc poison.cc util.cc +SRCS = attributes.cc distribute.cc input.cc iof5.cc output.cc poison.cc util.cc # Subdirectories containing source files SUBDIRS = diff --git a/CarpetDev/CarpetIOF5/src/output.cc b/CarpetDev/CarpetIOF5/src/output.cc index d376bbfcf..8a5ab5e53 100644 --- a/CarpetDev/CarpetIOF5/src/output.cc +++ b/CarpetDev/CarpetIOF5/src/output.cc @@ -41,10 +41,11 @@ namespace CarpetIOF5 { class output_iterator_t { // Can't be cGH const, since the mode loops change change its // entries - cGH* const cctkGH; + cGH *const cctkGH; vector<bool> const output_var; // whether to output this variable - bool const output_everything; + bool const output_past_timelevels; + bool const output_metadata; bool const is_multipatch; int group_type; // CCTK_GF or CCTK_ARRAY @@ -64,7 +65,7 @@ namespace CarpetIOF5 { rvect origin, delta; rvect lower, upper; - map_indices_t (cGH const* const cctkGH, int const gindex) + map_indices_t(cGH const *const cctkGH, int const gindex) { DECLARE_CCTK_ARGUMENTS; @@ -79,8 +80,8 @@ namespace CarpetIOF5 { } else { // grid array cGroupDynamicData dyndata; - int const ierr = CCTK_GroupDynamicData (cctkGH, gindex, &dyndata); - assert (not ierr); + int const ierr = CCTK_GroupDynamicData(cctkGH, gindex, &dyndata); + assert(not ierr); // HDF5 and F5 can't handle dim=0 dim = max(dyndata.dim, 1); for (int d=0; d<dyndata.dim; ++d) { @@ -103,11 +104,11 @@ namespace CarpetIOF5 { struct component_indices_t: map_indices_t { // elements >=dim remain undefined - ivect lbnd, lsh, lghosts, ughosts; + ivect lbnd, lssh, lghosts, ughosts; ivect imin, imax, ioff, ilen; rvect clower, cupper; - component_indices_t (cGH const* const cctkGH, int const gindex) + component_indices_t(cGH const *const cctkGH, int const gindex) : map_indices_t(cctkGH, gindex) { DECLARE_CCTK_ARGUMENTS; @@ -117,32 +118,37 @@ namespace CarpetIOF5 { // grid function for (int d=0; d<(::dim); ++d) { lbnd[d] = cctk_lbnd[d]; - lsh[d] = cctk_lsh[d]; + lssh[d] = CCTK_LSSH(0,d); lghosts[d] = cctk_bbox[2*d ] ? 0 : cctk_nghostzones[d]; ughosts[d] = cctk_bbox[2*d+1] ? 0 : cctk_nghostzones[d]; } } else { // grid array cGroupDynamicData dyndata; - int const ierr = CCTK_GroupDynamicData (cctkGH, gindex, &dyndata); - assert (not ierr); + int const ierr = CCTK_GroupDynamicData(cctkGH, gindex, &dyndata); + assert(not ierr); for (int d=0; d<dyndata.dim; ++d) { lbnd[d] = dyndata.lbnd[d]; - lsh[d] = dyndata.lsh[d]; +#ifdef CCTK_GROUPDYNAMICDATA_HAS_LSSH + lssh[d] = dyndata.lssh[CCTK_LSSH_IDX(0,d)]; +#else + lssh[d] = dyndata.lsh[d]; +#endif lghosts[d] = dyndata.bbox[2*d ] ? 0 : dyndata.nghostzones[d]; ughosts[d] = dyndata.bbox[2*d+1] ? 0 : dyndata.nghostzones[d]; } for (int d=dyndata.dim; d<(::dim); ++d) { lbnd[d] = 0; - lsh[d] = 1; + lssh[d] = 1; lghosts[d] = 0; ughosts[d] = 0; } } for (int d=0; d<(::dim); ++d) { imin[d] = 0; - imax[d] = lsh[d]; + imax[d] = lssh[d]; if (not output_ghost_points) { +#warning "TODO: Don't output ghosts on refinement boundaries; only output ghosts for inter-process boundaries" int const overlap = min(ughosts[d], minimum_component_overlap); imin[d] += lghosts[d]; imax[d] -= ughosts[d] - overlap; @@ -160,18 +166,20 @@ namespace CarpetIOF5 { public: - output_iterator_t (cGH* const cctkGH_, - vector<bool> const& output_var_, - bool const output_everything_) + output_iterator_t(cGH *const cctkGH_, + vector<bool> const& output_var_, + bool const output_past_timelevels_, + bool const output_metadata_) : cctkGH(cctkGH_), output_var(output_var_), - output_everything(output_everything_), + output_past_timelevels(output_past_timelevels_), + output_metadata(output_metadata_), is_multipatch (CCTK_IsFunctionAliased("MultiPatch_GetSystemSpecification")) { } - void iterate (hid_t const file) + void iterate(hid_t const file) { // Iterate over the variables in groups, first all grid // functions, then all non-GF groups @@ -182,28 +190,28 @@ namespace CarpetIOF5 { if (output_var.at(vindex)) { int const gindex = CCTK_GroupIndexFromVarI(vindex); if (CCTK_GroupTypeI(gindex) == CCTK_GF) { - vindices.push_back (vindex); + vindices.push_back(vindex); } } } if (not vindices.empty()) { - output_simulation (file); + output_simulation(file); } group_type = CCTK_ARRAY; for (group_index=0; group_index<CCTK_NumGroups(); ++group_index) { if (CCTK_GroupTypeI(group_index) != CCTK_GF) { vindices.clear(); - int const first_vindex = CCTK_FirstVarIndexI (group_index); - int const num_vars = CCTK_NumVarsInGroupI (group_index); + int const first_vindex = CCTK_FirstVarIndexI(group_index); + int const num_vars = CCTK_NumVarsInGroupI(group_index); for (int vindex=first_vindex; vindex<first_vindex+num_vars; ++vindex) { if (output_var.at(vindex)) { - vindices.push_back (vindex); + vindices.push_back(vindex); } } if (not vindices.empty()) { - output_simulation (file); + output_simulation(file); } } } @@ -213,39 +221,41 @@ namespace CarpetIOF5 { private: - void output_simulation (hid_t const file) + void output_simulation(hid_t const file) { DECLARE_CCTK_PARAMETERS; indent_t indent; - bool error_flag = false; gridname = generate_gridname(cctkGH); cout << indent << "simulation=" << gridname << "\n"; - assert (is_global_mode()); + assert(is_global_mode()); chartname = generate_chartname(cctkGH); + cout << indent << "chartname=" << chartname << "\n"; int const max_rl = group_type == CCTK_GF ? reflevels : 1; for (int rl=0; rl<max_rl; ++rl) { // Continue only if this level exists at this iteration - assert (maxtimereflevelfact % timereffacts.AT(rl) == 0); + assert(maxtimereflevelfact % timereffacts.AT(rl) == 0); int const do_every = group_type == CCTK_GF ? maxtimereflevelfact / timereffacts.AT(rl) : 1; if (cctkGH->cctk_iteration % do_every == 0) { +#warning "TODO: don't switch modes" ENTER_LEVEL_MODE(cctkGH, rl) { DECLARE_CCTK_ARGUMENTS; - assert (timelevel == 0); + assert(timelevel == 0); int const max_tl = - output_everything or output_all_timelevels ? + output_past_timelevels or output_all_timelevels ? (group_type == CCTK_GF ? timelevels : CCTK_NumTimeLevelsI(group_index)) : 1; for (timelevel=0; timelevel<max_tl; ++timelevel) { - cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel); + cctkGH->cctk_time = tt->get_time(mglevel, reflevel, timelevel); +#if 0 // Choose (arbitrarily) the root level as default // topology, for readers which don't understand AMR if (reflevel == 0) { @@ -253,43 +263,24 @@ namespace CarpetIOF5 { F5Path *const globalpath = F5Rcreate_vertex_refinement3D (file, cctk_time, gridname.c_str(), &v2h(reffact)[0], chartname.c_str()); - assert (globalpath); - FAILWARN (F5Rlink_default_vertex_topology (globalpath, - &v2h(reffact)[0])); + assert(globalpath); +#warning "TODO: Probably must not call this for cell-centred AMR; this probably makes the call to F5Rcreate_coordinate_topology below fail" + FAILWARN(F5Rlink_default_vertex_topology(globalpath, + &v2h(reffact)[0])); // Define iteration - FAILWARN (F5Rset_timestep (globalpath, cctk_iteration)); - - // Attach Cactus/Carpet metadata - if (timelevel == 0) { - // hid_t const metadata_group = globalpath->Grid_hid; - ostringstream pathname; - pathname << FIBER_CONTENT_GRIDS << "/" << gridname; - hid_t group; - group = H5Gopen (globalpath->ContentsGroup_hid, - pathname.str().c_str(), - H5P_DEFAULT); - if (group < 0) { - group = H5Gcreate (globalpath->ContentsGroup_hid, - pathname.str().c_str(), - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - } - assert (group >= 0); - write_metadata (cctkGH, group); - herr_t const herr = H5Gclose (group); - assert (not herr); - } + FAILWARN(F5Rset_timestep(globalpath, cctk_iteration)); // Close topology - F5close (globalpath); - - } // if reflevel==0 + F5close(globalpath); + } +#endif - output_reflevel (file); + output_reflevel(file); } // for timelevel timelevel = 0; - cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel); + cctkGH->cctk_time = tt->get_time(mglevel, reflevel, timelevel); } LEAVE_LEVEL_MODE; } // if do_every @@ -299,12 +290,13 @@ namespace CarpetIOF5 { - void output_reflevel (hid_t const file) + void output_reflevel(hid_t const file) { DECLARE_CCTK_ARGUMENTS; indent_t indent; + bool error_flag = false; - assert (is_level_mode()); + assert(is_level_mode()); ivect const reffact = spacereffacts.AT(reflevel); topologyname = generate_topologyname(cctkGH, group_index, reffact); @@ -314,26 +306,100 @@ namespace CarpetIOF5 { // Define grid hierarchy map_indices_t const mi(cctkGH, group_index); - int const indexdepth = vhh.at(0)->refcent == vertex_centered ? 0 : 1; - F5Path *const path = - F5Rcreate_coordinate_topology (file, &cctk_time, - gridname.c_str(), chartname.c_str(), - topologyname.c_str(), - indexdepth, - mi.dim, mi.dim, &v2h(reffact)[0]); - assert (path); - - BEGIN_LOCAL_MAP_LOOP (cctkGH, CCTK_GF) { - output_map (path); + F5Path *path = NULL; + F5Path *coordpath = NULL; + bool close_coordpath = true; + if (group_type == CCTK_GF) { + int const indexdepth = vhh.at(0)->refcent == vertex_centered ? 0 : 1; + if (indexdepth == 0) { + path = + F5Rcreate_coordinate_topology(file, &cctk_time, + gridname.c_str(), chartname.c_str(), + topologyname.c_str(), + 0, + mi.dim, mi.dim, &v2h(reffact)[0]); + } else { + char vertextopologyname[1000]; + TopologyName(vertextopologyname, sizeof vertextopologyname, + &v2h(reffact)[0], 0, dim); + coordpath = + F5Rcreate_coordinate_topology(file, &cctk_time, + gridname.c_str(), chartname.c_str(), + vertextopologyname, + 0, + mi.dim, mi.dim, &v2h(reffact)[0]); + path = + F5Rcreate_coordinate_topology(file, &cctk_time, + gridname.c_str(), + /*chartname.c_str(),*/ + vertextopologyname, + topologyname.c_str(), + indexdepth, + mi.dim, mi.dim, &v2h(reffact)[0]); +#warning "TODO: how should these two topologies be linked?" + // assert(mi.dim == 3); + // assert(all(reffact == 1)); + // path = + // F5Rcreate_hexaedrons_as_vertices_topology(file, &cctk_time, + // gridname.c_str()); + } + } else { + path = + F5Rcreate_coordinate_topology(file, &cctk_time, + gridname.c_str(), chartname.c_str(), + topologyname.c_str(), + 0, + mi.dim, mi.dim, &v2h(reffact)[0]); + } + if (!coordpath) { + coordpath = path; + close_coordpath = false; + } + assert(coordpath); + assert(path); + +#warning "TODO: Attach the number of I/O processes to the topology. WB sent email with suggestions." + + // Define default topology (once per grid) + if (group_type == CCTK_GF and reflevel == 0 and timelevel == 0) { + FAILWARN(F5Rlink_default_vertex_topology(path, &v2h(reffact)[0])); + + // Define iteration + FAILWARN(F5Rset_timestep(path, cctk_iteration)); + } + + // Attach Cactus/Carpet metadata (only once per output) + if (output_metadata) { + if (group_type == CCTK_GF and reflevel == 0 and timelevel == 0) { + // hid_t const metadata_group = path->Grid_hid; + ostringstream pathname; + pathname << FIBER_CONTENT_GRIDS << "/" << gridname; + hid_t group; + group = H5Gopen(path->ContentsGroup_hid, pathname.str().c_str(), + H5P_DEFAULT); + if (group < 0) { + group = H5Gcreate(path->ContentsGroup_hid, pathname.str().c_str(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + } + assert(group >= 0); + write_metadata(cctkGH, group); + herr_t const herr = H5Gclose(group); + assert(not herr); + } + } + + BEGIN_LOCAL_MAP_LOOP(cctkGH, group_type) { + output_map(path); } END_LOCAL_MAP_LOOP; - // Close topology - F5close (path); + // Close topologies + if (close_coordpath) F5close(coordpath); + F5close(path); } - void output_map (F5Path *const path) + void output_map(F5Path *const path) { DECLARE_CCTK_ARGUMENTS; indent_t indent; @@ -341,9 +407,9 @@ namespace CarpetIOF5 { cout << indent << "map=" << Carpet::map << "\n"; - assert (is_singlemap_mode()); + assert(is_singlemap_mode()); - if (not is_multipatch) { + if (group_type != CCTK_GF or not is_multipatch) { // Define level geometry map_indices_t const mi(cctkGH, group_index); F5_vec3_double_t const vlower = v2d(mi.lower); @@ -355,42 +421,42 @@ namespace CarpetIOF5 { charts.at(3) = F5B_standard_cartesian_chart3D(); } if (not charts.at(mi.dim)) { - assert (mi.dim != 0); - char const* coordnames[] = {"x", "y", "z"}; + assert(mi.dim != 0); + char const *coordnames[] = {"x", "y", "z"}; ostringstream chartnamebuf; chartnamebuf << "Cartesian " << mi.dim << "D"; charts.at(mi.dim) = F5B_new_global_float_chart(coordnames, mi.dim, chartnamebuf.str().c_str(), F5_FORTRAN_ORDER); - assert (charts.at(mi.dim)); + assert(charts.at(mi.dim)); } // hid_t const type = charts.at(mi.dim)->DoublePrecision.Point_hid_t; hid_t const type = path->myChart->DoublePrecision.Point_hid_t; - assert (type); - FAILWARN (F5Fwrite_linear (path, FIBER_HDF5_POSITIONS_STRING, - mi.dim, &v2h(mi.gsh)[0], - type, - &vlower, &vdelta)); + assert(type); + FAILWARN(F5Fwrite_linear(path, FIBER_HDF5_POSITIONS_STRING, + mi.dim, &v2h(mi.gsh)[0], + type, + &vlower, &vdelta)); #warning "TODO: path and chart don't match" - FAILWARN (F5Fset_range (path, &vlower, &vupper)); + FAILWARN(F5Fset_range(path, &vlower, &vupper)); } - BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) { - output_component (path); + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, group_type) { + output_component(path); } END_LOCAL_COMPONENT_LOOP; } - void output_component (F5Path *const path) + void output_component(F5Path *const path) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; indent_t indent; bool error_flag = false; - assert (is_local_mode()); + assert(is_local_mode()); fragmentname = generate_fragmentname(cctkGH, Carpet::map, component); cout << indent @@ -406,16 +472,16 @@ namespace CarpetIOF5 { // box was already defined above, but it provides the // individual components' bounding boxes.) component_indices_t const ci(cctkGH, group_index); - FAILWARN (F5Fwrite_linear_fraction (path, FIBER_HDF5_POSITIONS_STRING, - ci.dim, - &v2h(ci.gsh)[0], &v2h(ci.ilen)[0], - F5T_COORD3_DOUBLE, - &ci.clower, &ci.delta, - &v2h(ci.ioff)[0], - fragmentname.c_str())); + FAILWARN(F5Fwrite_linear_fraction(path, FIBER_HDF5_POSITIONS_STRING, + ci.dim, + &v2h(ci.gsh)[0], &v2h(ci.ilen)[0], + F5T_COORD3_DOUBLE, + &ci.clower, &ci.delta, + &v2h(ci.ioff)[0], + fragmentname.c_str())); } else { // Output coordinates - output_variable (path, CCTK_VarIndex("grid::x"), true); + output_variable(path, CCTK_VarIndex("grid::x"), true); } } @@ -423,14 +489,14 @@ namespace CarpetIOF5 { for (vector<int>::const_iterator vi = vindices.begin(); vi != vindices.end(); ++vi) { - output_variable (path, *vi); + output_variable(path, *vi); } } - void output_variable (F5Path *const path, int const var, - bool const write_positions = false) + void output_variable(F5Path *const path, int const var, + bool const write_positions = false) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; @@ -439,7 +505,7 @@ namespace CarpetIOF5 { int ierr; - assert (var >= 0); + assert(var >= 0); if (write_positions) { cout << indent << "positions\n"; } else { @@ -447,18 +513,21 @@ namespace CarpetIOF5 { char *const groupname = CCTK_GroupNameFromVarI(var); cout << indent << "variable=" << fullname << " (group=" << groupname << ")\n"; - free (groupname); - free (fullname); + free(groupname); + free(fullname); } int const group = CCTK_GroupIndexFromVarI(var); int const v0 = CCTK_FirstVarIndexI(group); + // Don't output groups without storage + if (not CCTK_QueryGroupStorageI(cctkGH, group)) return; + cGroup groupdata; - ierr = CCTK_GroupData (group, &groupdata); - assert (not ierr); + ierr = CCTK_GroupData(group, &groupdata); + assert(not ierr); - assert ((groupdata.grouptype == CCTK_GF) == (group_type == CCTK_GF)); + assert((groupdata.grouptype == CCTK_GF) == (group_type == CCTK_GF)); // Output distrib=constant variables only on process 0 switch (groupdata.disttype) { @@ -469,19 +538,20 @@ namespace CarpetIOF5 { // do nothing break; default: - assert (0); + assert(0); } - assert (groupdata.stagtype == 0); + assert(groupdata.stagtype == 0); #warning "TODO: Do not output symmetry zones (unless requested by the user)" #warning "TODO: Do not output buffer zones (is that easily possible?)" - int const will_cover_complete_domain = not is_multipatch and reflevel==0; + int const will_cover_complete_domain = + (group_type != CCTK_GF or not is_multipatch) and reflevel==0; cGroupDynamicData dyndata; - ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata); - assert (not ierr); + ierr = CCTK_GroupDynamicData(cctkGH, group, &dyndata); + assert(not ierr); // Only output active timelevels if (timelevel >= dyndata.activetimelevels) return; @@ -490,8 +560,8 @@ namespace CarpetIOF5 { tensortype_t tensortype = tt_error; - int const coordinates_group = CCTK_GroupIndex ("grid::coordinates"); - assert (coordinates_group >= 0); + int const coordinates_group = CCTK_GroupIndex("grid::coordinates"); + assert(coordinates_group >= 0); if (group == coordinates_group) { @@ -501,7 +571,7 @@ namespace CarpetIOF5 { } else if (var == v0+dim) { tensortype = tt_scalar; } else { - assert (0); + assert(0); } } else { @@ -541,14 +611,14 @@ namespace CarpetIOF5 { case tt_scalar: { // Scalar field - assert (not write_positions); + assert(not write_positions); switch (groupdata.vartype) { case CCTK_VARIABLE_INT: type = H5T_NATIVE_INT; break; case CCTK_VARIABLE_REAL: type = H5T_NATIVE_DOUBLE; break; default: assert(0); } num_comps = 1; - name = generate_fieldname (cctkGH, var, tensortype); + name = generate_fieldname(cctkGH, var, tensortype); break; } @@ -564,17 +634,17 @@ namespace CarpetIOF5 { name = write_positions ? FIBER_HDF5_POSITIONS_STRING : - generate_fieldname (cctkGH, var, tensortype); + generate_fieldname(cctkGH, var, tensortype); if (write_positions) { htri_t const exists = - H5Lexists (path->Representation_hid, name.c_str(), H5P_DEFAULT); - assert (exists >= 0); + H5Lexists(path->Representation_hid, name.c_str(), H5P_DEFAULT); + assert(exists >= 0); if (exists) { string const fragmentpath = name + "/" + fragmentname; htri_t const exists2 = - H5Lexists (path->Representation_hid, fragmentpath.c_str(), - H5P_DEFAULT); - assert (exists2 >= 0); + H5Lexists(path->Representation_hid, fragmentpath.c_str(), + H5P_DEFAULT); + assert(exists2 >= 0); if (exists2) { // Write positions only once indent_t indent2; @@ -588,43 +658,46 @@ namespace CarpetIOF5 { case tt_symtensor: { // Symmetric tensor field - assert (not write_positions); + assert(not write_positions); switch (groupdata.vartype) { case CCTK_VARIABLE_REAL: type = F5T_METRIC33_DOUBLE; break; default: assert(0); } num_comps = dim*(dim+1)/2; - name = generate_fieldname (cctkGH, var, tensortype); + name = generate_fieldname(cctkGH, var, tensortype); break; } case tt_tensor: { // Non-symmetric tensor field - assert (not write_positions); + assert(not write_positions); switch (groupdata.vartype) { case CCTK_VARIABLE_REAL: type = F5T_BIVEC3_DOUBLE; break; default: assert(0); } num_comps = dim*dim; - name = generate_fieldname (cctkGH, var, tensortype); + name = generate_fieldname(cctkGH, var, tensortype); break; } default: - assert (0); + assert(0); } cout << indent << "fieldname=" << name << "\n"; // Write data - assert (type >= 0); - assert (num_comps > 0); - assert (name.size() > 0); + assert(type >= 0); + assert(num_comps > 0); + assert(name.size() > 0); + + // Ensure that the data types match + int const vartype = CCTK_VarTypeI(var); + assert(num_comps * CCTK_VarTypeSize(vartype) == (int)H5Tget_size(type)); component_indices_t const ci(cctkGH, group_index); - void const* data[num_comps]; - int const vartype = CCTK_VarTypeI(var); + void const *data[num_comps]; switch (vartype) { case CCTK_VARIABLE_INT: { @@ -637,8 +710,8 @@ namespace CarpetIOF5 { for (int j=0; j<ci.ilen[1]; ++j) { for (int i=0; i<ci.ilen[0]; ++i) { int const isrc = - (ci.imin[0]+i + ci.lsh[0] * - (ci.imin[1]+j + ci.lsh[1] * + (ci.imin[0]+i + ci.lssh[0] * + (ci.imin[1]+j + ci.lssh[1] * (ci.imin[2]+k))); int const idst = i + ci.ilen[0] * (j + ci.ilen[1] * k); rdata[idst] = varptr[isrc]; @@ -660,8 +733,8 @@ namespace CarpetIOF5 { for (int j=0; j<ci.ilen[1]; ++j) { for (int i=0; i<ci.ilen[0]; ++i) { int const isrc = - (ci.imin[0]+i + ci.lsh[0] * - (ci.imin[1]+j + ci.lsh[1] * + (ci.imin[0]+i + ci.lssh[0] * + (ci.imin[1]+j + ci.lssh[1] * (ci.imin[2]+k))); int const idst = i + ci.ilen[0] * (j + ci.ilen[1] * k); rdata[idst] = varptr[isrc]; @@ -678,17 +751,17 @@ namespace CarpetIOF5 { } // Dataset properties - hid_t const prop = H5Pcreate (H5P_DATASET_CREATE); - assert (prop >= 0); + hid_t const prop = H5Pcreate(H5P_DATASET_CREATE); + assert(prop >= 0); // Enable compression if requested if (compression_level >= 0) { - FAILWARN (H5Pset_chunk (prop, ci.dim, &v2h(ci.ilen).reverse()[0])); - FAILWARN (H5Pset_deflate (prop, compression_level)); + FAILWARN(H5Pset_chunk(prop, ci.dim, &v2h(ci.ilen).reverse()[0])); + FAILWARN(H5Pset_deflate(prop, compression_level)); } // Enable checksums if requested if (use_checksums) { - FAILWARN (H5Pset_chunk (prop, ci.dim, &v2h(ci.ilen).reverse()[0])); - FAILWARN (H5Pset_fletcher32 (prop)); + FAILWARN(H5Pset_chunk(prop, ci.dim, &v2h(ci.ilen).reverse()[0])); + FAILWARN(H5Pset_fletcher32(prop)); } if (num_comps == 1 and not separate_single_component_tensors) { @@ -697,27 +770,39 @@ namespace CarpetIOF5 { // instead) // TODO: Extent F5 API to allow writing non-fragmented datasets FAILWARN - (F5Fwrite_fraction (path, name.c_str(), + (F5Fwrite_fraction(path, name.c_str(), + ci.dim, + is_multipatch ? NULL : &v2h(ci.gsh)[0], + &v2h(ci.ilen)[0], + type, type, data[0], + &v2h(ci.ioff)[0], + &v2h(ci.lghosts)[0], &v2h(ci.ughosts)[0], + fragmentname.c_str(), prop)); + } else { + int const full_coverage = + will_cover_complete_domain and not fragment_contiguous_components; + // F5ls does not seem to support what F5 generates in this + // case. It seems that F5 does not generate a separated object + // (it does not generate a group for all datasets), but + // declares that it does create one (it calls it + // FragmentedSeparatedCompound). F5ls then gets confused + // because it cannot open this group. + if (not full_coverage and num_comps == 1) { + CCTK_WARN(CCTK_WARN_ALERT, + "Outputting scalars in a fragmented, separated way. " + "This does not seem to be supported by F5ls " + "(or is implemented wrong in the F5 library)."); + } + FAILWARN + (F5FSwrite_fraction(path, name.c_str(), ci.dim, is_multipatch ? NULL : &v2h(ci.gsh)[0], &v2h(ci.ilen)[0], - type, type, data[0], + type, type, data, &v2h(ci.ioff)[0], &v2h(ci.lghosts)[0], &v2h(ci.ughosts)[0], - fragmentname.c_str(), prop)); - } else { - int const full_coverage = - will_cover_complete_domain and not fragment_contiguous_components; - FAILWARN - (F5FSwrite_fraction (path, name.c_str(), - ci.dim, - is_multipatch ? NULL : &v2h(ci.gsh)[0], - &v2h(ci.ilen)[0], - type, type, data, - &v2h(ci.ioff)[0], - &v2h(ci.lghosts)[0], &v2h(ci.ughosts)[0], - fragmentname.c_str(), prop, - full_coverage)); + fragmentname.c_str(), prop, + full_coverage)); } for (int d=0; d<num_comps; ++d) { @@ -729,7 +814,7 @@ namespace CarpetIOF5 { data[d] = NULL; } - FAILWARN (H5Pclose (prop)); + FAILWARN(H5Pclose(prop)); } @@ -737,50 +822,47 @@ namespace CarpetIOF5 { - void write_metadata (cGH const* const cctkGH, hid_t const file) + void write_metadata(cGH const *const cctkGH, hid_t const file) { DECLARE_CCTK_PARAMETERS; - assert (cctkGH); - assert (file >= 0); + assert(cctkGH); + assert(file >= 0); herr_t herr; - CCTK_INFO ("Writing simulation metadata..."); - // NOTE: We could write the metadata into only one of the process - // files (to save space), or write it into a separate metadata - // file + CCTK_INFO("Writing simulation metadata..."); // Create metadata only once - // TODO: update metadata at some point - htri_t const exists = H5Lexists (file, metadata_group, H5P_DEFAULT); - assert (exists >= 0); + // TODO: instead, overwrite the metadata + htri_t const exists = H5Lexists(file, metadata_group, H5P_DEFAULT); + assert(exists >= 0); if (exists) return; // Create a group to hold all metadata hid_t const group = - H5Gcreate (file, metadata_group, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - assert (group >= 0); + H5Gcreate(file, metadata_group, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(group >= 0); // General information - WriteAttribute (group, "Cactus version", CCTK_FullVersion()); + WriteAttribute(group, "Cactus version", CCTK_FullVersion()); // Unique identifiers - if (CCTK_IsFunctionAliased ("UniqueConfigID")) { + if (CCTK_IsFunctionAliased("UniqueConfigID")) { WriteAttribute (group, "config id", (char const*) UniqueConfigID(cctkGH)); } - if (CCTK_IsFunctionAliased ("UniqueBuildID")) { + if (CCTK_IsFunctionAliased("UniqueBuildID")) { WriteAttribute (group, "build id", (char const*) UniqueBuildID(cctkGH)); } - if (CCTK_IsFunctionAliased ("UniqueSimulationID")) { + if (CCTK_IsFunctionAliased("UniqueSimulationID")) { WriteAttribute (group, "simulation id", (char const*) UniqueSimulationID(cctkGH)); } - if (CCTK_IsFunctionAliased ("UniqueRunID")) { + if (CCTK_IsFunctionAliased("UniqueRunID")) { WriteAttribute (group, "run id", (char const*) UniqueRunID(cctkGH)); } @@ -791,35 +873,37 @@ namespace CarpetIOF5 { // Number of I/O processes (i.e. the number of output files) int const nprocs = CCTK_nProcs(cctkGH); int const nioprocs = nprocs; - WriteAttribute (group, "nioprocs", nioprocs); + WriteAttribute(group, "nioprocs", nioprocs); #endif // Write parameters into a separate dataset (they may be too large // for an attribute) int const get_all = 1; - char *const parameters = IOUtil_GetAllParameters (cctkGH, get_all); - assert (parameters); - WriteLargeAttribute (group, all_parameters, parameters); - free (parameters); + char *const parameters = IOUtil_GetAllParameters(cctkGH, get_all); + assert(parameters); + WriteLargeAttribute(group, all_parameters, parameters); + free(parameters); // Grid structure - string const gs = serialise_grid_structure (cctkGH); - WriteLargeAttribute (group, grid_structure, gs.c_str()); + string const gs = serialise_grid_structure(cctkGH); + WriteLargeAttribute(group, grid_structure, gs.c_str()); - herr = H5Gclose (group); - assert (not herr); + herr = H5Gclose(group); + assert(not herr); } - void output (cGH const* const cctkGH, - hid_t const file, - vector<bool> const& output_var, - bool const output_everything) + void output(cGH const *const cctkGH, + hid_t const file, + vector<bool> const& output_var, + bool const output_past_timelevels, + bool const output_metadata) { output_iterator_t - iterator(const_cast<cGH*>(cctkGH), output_var, output_everything); - iterator.iterate (file); + iterator(const_cast<cGH*>(cctkGH), output_var, + output_past_timelevels, output_metadata); + iterator.iterate(file); } } // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/poison.cc b/CarpetDev/CarpetIOF5/src/poison.cc index 8dcc9cbdb..f9292def5 100644 --- a/CarpetDev/CarpetIOF5/src/poison.cc +++ b/CarpetDev/CarpetIOF5/src/poison.cc @@ -17,13 +17,13 @@ namespace CarpetIOF5 { extern "C" - void F5_Poison (CCTK_ARGUMENTS) + void F5_Poison(CCTK_ARGUMENTS) { DECLARE_CCTK_PARAMETERS; - assert (is_global_mode()); - CCTK_VInfo (CCTK_THORNSTRING, - "F5_Poison: iteration=%d", cctkGH->cctk_iteration); + assert(is_global_mode()); + CCTK_VInfo(CCTK_THORNSTRING, + "F5_Poison: iteration=%d", cctkGH->cctk_iteration); BEGIN_REFLEVEL_LOOP(cctkGH) { BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) { @@ -45,13 +45,13 @@ namespace CarpetIOF5 { extern "C" - void F5_Check (CCTK_ARGUMENTS) + void F5_Check(CCTK_ARGUMENTS) { DECLARE_CCTK_PARAMETERS; - assert (is_global_mode()); - CCTK_VInfo (CCTK_THORNSTRING, - "F5_Check: iteration=%d", cctkGH->cctk_iteration); + assert(is_global_mode()); + CCTK_VInfo(CCTK_THORNSTRING, + "F5_Check: iteration=%d", cctkGH->cctk_iteration); BEGIN_REFLEVEL_LOOP(cctkGH) { BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) { @@ -60,10 +60,10 @@ namespace CarpetIOF5 { #pragma omp parallel CCTK_LOOP3_ALL(F5_Check, cctkGH, i,j,k) { int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k); - assert (not isnan(x[ind3d])); - assert (not isnan(y[ind3d])); - assert (not isnan(z[ind3d])); - assert (not isnan(r[ind3d])); + assert(not isnan(x[ind3d])); + assert(not isnan(y[ind3d])); + assert(not isnan(z[ind3d])); + assert(not isnan(r[ind3d])); } CCTK_ENDLOOP3_ALL(F5_Check); } END_LOCAL_COMPONENT_LOOP; } END_LOCAL_MAP_LOOP; diff --git a/CarpetDev/CarpetIOF5/src/util.cc b/CarpetDev/CarpetIOF5/src/util.cc index ce3c83f5e..3d924e775 100644 --- a/CarpetDev/CarpetIOF5/src/util.cc +++ b/CarpetDev/CarpetIOF5/src/util.cc @@ -67,61 +67,61 @@ namespace CarpetIOF5 { // Generate a good file name ("alias") for a variable string - generate_basename (cGH const *const cctkGH, - int const variable) + generate_basename(cGH const *const cctkGH, + int const variable) { DECLARE_CCTK_PARAMETERS; ostringstream filename_buf; - if (CCTK_EQUALS (file_content, "variable")) { - assert (variable >= 0); + if (CCTK_EQUALS(file_content, "variable")) { + assert(variable >= 0); char *const varname = CCTK_FullName(variable); - assert (varname); + assert(varname); for (char *p = varname; *p; ++p) { *p = tolower(*p); } filename_buf << varname; - free (varname); + free(varname); } - else if (CCTK_EQUALS (file_content, "group")) { - assert (variable >= 0); + else if (CCTK_EQUALS(file_content, "group")) { + assert(variable >= 0); char *const groupname = CCTK_GroupNameFromVarI(variable); - assert (groupname); + assert(groupname); for (char *p = groupname; *p; ++p) { *p = tolower(*p); } filename_buf << groupname; - free (groupname); + free(groupname); } - else if (CCTK_EQUALS (file_content, "thorn")) { - assert (variable >= 0); + else if (CCTK_EQUALS(file_content, "thorn")) { + assert(variable >= 0); char const *const impname = CCTK_ImpFromVarI(variable); char *const thornname = strdup(impname); - assert (thornname); + assert(thornname); char *const colon = strchr(thornname, ':'); - assert (colon); + assert(colon); *colon = '\0'; for (char *p = thornname; *p; ++p) { *p = tolower(*p); } filename_buf << thornname; - free (thornname); + free(thornname); } - else if (CCTK_EQUALS (file_content, "everything")) { + else if (CCTK_EQUALS(file_content, "everything")) { if (CCTK_EQUALS(out_filename, "")) { // Obtain the parameter file name char buf[10000]; - int const ilen = CCTK_ParameterFilename (sizeof buf, buf); - assert (ilen < int(sizeof buf) - 1); - char* parfilename = buf; + int const ilen = CCTK_ParameterFilename(sizeof buf, buf); + assert(ilen < int(sizeof buf) - 1); + char *parfilename = buf; // Remove directory prefix, if any - char* const slash = strrchr(parfilename, '/'); + char *const slash = strrchr(parfilename, '/'); if (slash) { parfilename = &slash[1]; } // Remove suffix, if it is there - char* const suffix = strrchr(parfilename, '.'); + char *const suffix = strrchr(parfilename, '.'); if (suffix and strcmp(suffix, ".par")==0) { suffix[0] = '\0'; } @@ -131,14 +131,14 @@ namespace CarpetIOF5 { } } else { - assert (0); + assert(0); } if (out_timesteps_per_file > 0) { int const iteration = (cctkGH->cctk_iteration / out_timesteps_per_file * out_timesteps_per_file); filename_buf << ".it" - << setw (iteration_digits) << setfill ('0') << iteration; + << setw(iteration_digits) << setfill('0') << iteration; } return filename_buf.str(); @@ -148,11 +148,12 @@ namespace CarpetIOF5 { // Create the final file name on a particular processor string - create_filename (cGH const *const cctkGH, - string const basename, - int const proc, - io_dir_t const io_dir, - bool const create_directories) + create_filename(cGH const *const cctkGH, + string const basename, + int const iteration, + int const proc, + io_dir_t const io_dir, + bool const create_directories) { DECLARE_CCTK_PARAMETERS; @@ -185,26 +186,26 @@ namespace CarpetIOF5 { if (create_subdirs) { { ostringstream buf; - buf << path + "/" + buf << path << "/" << basename - << ".p" << setw (max (0, processor_digits - 4)) << setfill ('0') + << ".p" << setw(max(0, processor_digits - 4)) << setfill('0') << proc / 10000 << "nnnn/"; path = buf.str(); if (create_directories) { - check (CCTK_CreateDirectory (mode, path.c_str()) >= 0); + check(CCTK_CreateDirectory(mode, path.c_str()) >= 0); } } { ostringstream buf; - buf << path + "/" + buf << path << "/" << basename - << ".p" << setw (max (0, processor_digits - 2)) << setfill ('0') + << ".p" << setw(max(0, processor_digits - 2)) << setfill('0') << proc / 100 << "nn/"; path = buf.str(); if (create_directories) { - check (CCTK_CreateDirectory (mode, path.c_str()) >= 0); + check(CCTK_CreateDirectory(mode, path.c_str()) >= 0); } } } @@ -212,38 +213,41 @@ namespace CarpetIOF5 { ostringstream buf; buf << path << basename - << ".p" << setw (processor_digits) << setfill ('0') + << ".p" << setw(processor_digits) << setfill('0') << proc << "/"; path = buf.str(); if (create_directories) { - check (CCTK_CreateDirectory (mode, path.c_str()) >= 0); + check(CCTK_CreateDirectory(mode, path.c_str()) >= 0); } } ostringstream buf; - buf << path + "/" - << basename - << ".p" << setw (processor_digits) << setfill ('0') << proc - << out_extension; - return buf.str(); + buf << path << "/" << basename; + if (io_dir==io_dir_recover or io_dir==io_dir_checkpoint) { + buf << ".i" << setw(iteration_digits) << setfill('0') << iteration; + } + buf << ".p" << setw(processor_digits) << setfill('0') << proc; + buf << out_extension; + path = buf.str(); + return path; } // Generate a good grid name (simulation name) string - generate_gridname (cGH const *const cctkGH) + generate_gridname(cGH const *const cctkGH) { #if 0 char const *const gridname = (char const*) UniqueSimulationID(cctkGH); - assert (gridname); + assert(gridname); return gridname; #endif // Use the parameter file name, since the simulation ID is too // long and doesn't look nice char parfilename[10000]; - CCTK_ParameterFilename (sizeof parfilename, parfilename); + CCTK_ParameterFilename(sizeof parfilename, parfilename); char const *const suffix = ".par"; if (strlen(parfilename) >= strlen(suffix)) { char *const p = parfilename + strlen(parfilename) - strlen(suffix); @@ -259,18 +263,19 @@ namespace CarpetIOF5 { // Generate a good topology name (refinement level name) string - generate_topologyname (cGH const *const cctkGH, - int const gi, - ivect const& reffact) + generate_topologyname(cGH const *const cctkGH, + int const gi, + ivect const& reffact) { ostringstream buf; if (gi == -1) { // grid function - buf << "VertexLevel_"; - for (int d=0; d<dim; ++d) { - if (d>0) buf << "x"; - buf << reffact[d]; - } + // bit d of centering indicates whether direction d is centered + int const centering = + vhh.at(0)->refcent == vertex_centered ? 0 : (1<<dim)-1; + char name[1000]; + TopologyName(name, sizeof name, &v2h(reffact)[0], centering, dim); + buf << name; } else { // grid array char *const groupname = CCTK_GroupName(gi); @@ -278,14 +283,14 @@ namespace CarpetIOF5 { *p = tolower(*p); } buf << "Group_" << groupname; - free (groupname); + free(groupname); } return buf.str(); } // Generate a good chart name (tensor basis name) string - generate_chartname (cGH const *const cctkGH) + generate_chartname(cGH const *const cctkGH) { // Assume a global tensor basis, where all maps use the same chart return FIBER_HDF5_DEFAULT_CHART; @@ -303,7 +308,7 @@ namespace CarpetIOF5 { // (We assume that fragment names need to be unique only within a // topology, not across topologies.) string - generate_fragmentname (cGH const *const cctkGH, int const m, int const c) + generate_fragmentname(cGH const *const cctkGH, int const m, int const c) { ostringstream buf; buf << "Map_" << m << "_Component_" << c; @@ -311,21 +316,21 @@ namespace CarpetIOF5 { } void - interpret_fragmentname (cGH const *const cctkGH, - char const *const fragmentname, - int& m, int& c) + interpret_fragmentname(cGH const *const cctkGH, + char const *const fragmentname, + int& m, int& c) { m = -1; c = -1; - sscanf (fragmentname, "Map_%d_Component_%d", &m, &c); - assert (m>=0 and m<Carpet::maps); - assert (c>=0); + sscanf(fragmentname, "Map_%d_Component_%d", &m, &c); + assert(m>=0 and m<Carpet::maps); + assert(c>=0); } // Generate a good field name (group or variable name) string - generate_fieldname (cGH const *const cctkGH, - int const vi, tensortype_t const tt) + generate_fieldname(cGH const *const cctkGH, + int const vi, tensortype_t const tt) { int const gi = CCTK_GroupIndexFromVarI(vi); int const numvars = CCTK_NumVarsInGroupI(gi); @@ -335,17 +340,17 @@ namespace CarpetIOF5 { if (tt == tt_scalar and numvars >1) { char *const fullname = CCTK_FullName(vi); name = fullname; - free (fullname); + free(fullname); } else { char *const groupname = CCTK_GroupName(gi); name = groupname; - free (groupname); + free(groupname); } - transform (name.begin(), name.end(), name.begin(), ::tolower); + transform(name.begin(), name.end(), name.begin(), ::tolower); string const sep = "::"; size_t const pos = name.find(sep); - assert (pos != string::npos); - name.replace (pos, sep.size(), "."); + assert(pos != string::npos); + name.replace(pos, sep.size(), "."); return name; } @@ -359,7 +364,7 @@ namespace CarpetIOF5 { vi = -1; return; } - fieldname.replace (pos, sep.size(), "::"); + fieldname.replace(pos, sep.size(), "::"); vi = CCTK_VarIndex(fieldname.c_str()); if (vi < 0) { @@ -370,7 +375,7 @@ namespace CarpetIOF5 { return; } vi = CCTK_FirstVarIndexI(gi); - assert (vi>=0); + assert(vi>=0); } } @@ -379,7 +384,7 @@ namespace CarpetIOF5 { char const *const grid_structure = "Grid Structure v5"; string - serialise_grid_structure (cGH const *const cctkGH) + serialise_grid_structure(cGH const *const cctkGH) { ostringstream buf; buf << setprecision(17); |