From 5fe6689c15dc6e5485dff1ecee016ff511f96cbe Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Tue, 24 Aug 2010 12:10:36 -0400 Subject: CarpetIOF5: New thorn --- CarpetDev/CarpetIOF5/README | 9 + CarpetDev/CarpetIOF5/configuration.ccl | 3 + CarpetDev/CarpetIOF5/doc/documentation.tex | 144 ++++++ CarpetDev/CarpetIOF5/interface.ccl | 46 ++ CarpetDev/CarpetIOF5/par/iof5-multipatch.par | 45 ++ CarpetDev/CarpetIOF5/par/iof5-refined.par | 48 ++ CarpetDev/CarpetIOF5/par/iof5-uniform.par | 39 ++ CarpetDev/CarpetIOF5/param.ccl | 66 +++ CarpetDev/CarpetIOF5/schedule.ccl | 27 ++ CarpetDev/CarpetIOF5/src/attributes.cc | 184 ++++++++ CarpetDev/CarpetIOF5/src/input.cc | 563 +++++++++++++++++++++++ CarpetDev/CarpetIOF5/src/iof5.cc.old | 643 +++++++++++++++++++++++++++ CarpetDev/CarpetIOF5/src/iof5.hh | 214 +++++++++ CarpetDev/CarpetIOF5/src/make.code.defn | 7 + CarpetDev/CarpetIOF5/src/output.cc | 573 ++++++++++++++++++++++++ CarpetDev/CarpetIOF5/src/poison.cc | 77 ++++ CarpetDev/CarpetIOF5/src/util.cc | 282 ++++++++++++ 17 files changed, 2970 insertions(+) create mode 100644 CarpetDev/CarpetIOF5/README create mode 100644 CarpetDev/CarpetIOF5/configuration.ccl create mode 100644 CarpetDev/CarpetIOF5/doc/documentation.tex create mode 100644 CarpetDev/CarpetIOF5/interface.ccl create mode 100644 CarpetDev/CarpetIOF5/par/iof5-multipatch.par create mode 100644 CarpetDev/CarpetIOF5/par/iof5-refined.par create mode 100644 CarpetDev/CarpetIOF5/par/iof5-uniform.par create mode 100644 CarpetDev/CarpetIOF5/param.ccl create mode 100644 CarpetDev/CarpetIOF5/schedule.ccl create mode 100644 CarpetDev/CarpetIOF5/src/attributes.cc create mode 100644 CarpetDev/CarpetIOF5/src/input.cc create mode 100644 CarpetDev/CarpetIOF5/src/iof5.cc.old create mode 100644 CarpetDev/CarpetIOF5/src/iof5.hh create mode 100644 CarpetDev/CarpetIOF5/src/make.code.defn create mode 100644 CarpetDev/CarpetIOF5/src/output.cc create mode 100644 CarpetDev/CarpetIOF5/src/poison.cc create mode 100644 CarpetDev/CarpetIOF5/src/util.cc (limited to 'CarpetDev') diff --git a/CarpetDev/CarpetIOF5/README b/CarpetDev/CarpetIOF5/README new file mode 100644 index 000000000..d17fddf2a --- /dev/null +++ b/CarpetDev/CarpetIOF5/README @@ -0,0 +1,9 @@ +Cactus Code Thorn CarpetIOF5 +Author(s) : Erik Schnetter +Maintainer(s): Erik Schnetter +Licence : GPL +-------------------------------------------------------------------------- + +1. Purpose + +Efficient, binary I/O based on the F5 library diff --git a/CarpetDev/CarpetIOF5/configuration.ccl b/CarpetDev/CarpetIOF5/configuration.ccl new file mode 100644 index 000000000..155a0c55b --- /dev/null +++ b/CarpetDev/CarpetIOF5/configuration.ccl @@ -0,0 +1,3 @@ +# Configuration definition for thorn CarpetIOF5 + +REQUIRES Carpet CarpetLib F5 HDF5 LoopControl diff --git a/CarpetDev/CarpetIOF5/doc/documentation.tex b/CarpetDev/CarpetIOF5/doc/documentation.tex new file mode 100644 index 000000000..e73a6e7ed --- /dev/null +++ b/CarpetDev/CarpetIOF5/doc/documentation.tex @@ -0,0 +1,144 @@ +% *======================================================================* +% Cactus Thorn template for ThornGuide documentation +% Author: Ian Kelley +% Date: Sun Jun 02, 2002 +% $Header$ +% +% Thorn documentation in the latex file doc/documentation.tex +% will be included in ThornGuides built with the Cactus make system. +% The scripts employed by the make system automatically include +% pages about variables, parameters and scheduling parsed from the +% relevant thorn CCL files. +% +% This template contains guidelines which help to assure that your +% documentation will be correctly added to ThornGuides. More +% information is available in the Cactus UsersGuide. +% +% Guidelines: +% - Do not change anything before the line +% % START CACTUS THORNGUIDE", +% except for filling in the title, author, date, etc. fields. +% - Each of these fields should only be on ONE line. +% - Author names should be separated with a \\ or a comma. +% - You can define your own macros, but they must appear after +% the START CACTUS THORNGUIDE line, and must not redefine standard +% latex commands. +% - To avoid name clashes with other thorns, 'labels', 'citations', +% 'references', and 'image' names should conform to the following +% convention: +% ARRANGEMENT_THORN_LABEL +% For example, an image wave.eps in the arrangement CactusWave and +% thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps +% - Graphics should only be included using the graphicx package. +% More specifically, with the "\includegraphics" command. Do +% not specify any graphic file extensions in your .tex file. This +% will allow us to create a PDF version of the ThornGuide +% via pdflatex. +% - References should be included with the latex "\bibitem" command. +% - Use \begin{abstract}...\end{abstract} instead of \abstract{...} +% - Do not use \appendix, instead include any appendices you need as +% standard sections. +% - For the benefit of our Perl scripts, and for future extensions, +% please use simple latex. +% +% *======================================================================* +% +% Example of including a graphic image: +% \begin{figure}[ht] +% \begin{center} +% \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure} +% \end{center} +% \caption{Illustration of this and that} +% \label{MyArrangement_MyThorn_MyLabel} +% \end{figure} +% +% Example of using a label: +% \label{MyArrangement_MyThorn_MyLabel} +% +% Example of a citation: +% \cite{MyArrangement_MyThorn_Author99} +% +% Example of including a reference +% \bibitem{MyArrangement_MyThorn_Author99} +% {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), +% 1--16. {\tt http://www.nowhere.com/}} +% +% *======================================================================* + +% If you are using CVS use this line to give version information +% $Header$ + +\documentclass{article} + +% Use the Cactus ThornGuide style file +% (Automatically used from Cactus distribution, if you have a +% thorn without the Cactus Flesh download this from the Cactus +% homepage at www.cactuscode.org) +\usepackage{../../../../doc/latex/cactus} + +\begin{document} + +% The author of the documentation +\author{Erik Schnetter \textless schnetter@cct.lsu.edu\textgreater} + +% The title of the document (not necessarily the name of the Thorn) +\title{CarpetIOF5} + +% the date your document was last changed, if your document is in CVS, +% please use: +% \date{$ $Date: 2004-01-07 14:12:39 -0600 (Wed, 07 Jan 2004) $ $} +\date{July 16 2010} + +\maketitle + +% Do not delete next line +% START CACTUS THORNGUIDE + +% Add all definitions used in this documentation here +% \def\mydef etc + +% Add an abstract for this thorn's documentation +\begin{abstract} + +\end{abstract} + +% The following sections are suggestive only. +% Remove them or add your own. + +\section{Introduction} + +\section{Physical System} + +\section{Numerical Implementation} + +\section{Using This Thorn} + +\subsection{Obtaining This Thorn} + +\subsection{Basic Usage} + +\subsection{Special Behaviour} + +\subsection{Interaction With Other Thorns} + +\subsection{Examples} + +\subsection{Support and Feedback} + +\section{History} + +\subsection{Thorn Source Code} + +\subsection{Thorn Documentation} + +\subsection{Acknowledgements} + + +\begin{thebibliography}{9} + +\end{thebibliography} + +% Do not delete next line +% END CACTUS THORNGUIDE + +\end{document} diff --git a/CarpetDev/CarpetIOF5/interface.ccl b/CarpetDev/CarpetIOF5/interface.ccl new file mode 100644 index 000000000..4c09d8b12 --- /dev/null +++ b/CarpetDev/CarpetIOF5/interface.ccl @@ -0,0 +1,46 @@ +# Interface definition for thorn CarpetIOF5 + +IMPLEMENTS: CarpetIOF5 + +INHERITS: grid + +USES INCLUDE: carpet.hh + + + +# Return a pointer to an unmodifiable C string +# which contains a unique ID for this configuration +CCTK_POINTER_TO_CONST \ +FUNCTION UniqueConfigID (CCTK_POINTER_TO_CONST IN cctkGH) +USES FUNCTION UniqueConfigID + +# Return a pointer to an unmodifiable C string +# which contains a unique ID for this build +CCTK_POINTER_TO_CONST \ +FUNCTION UniqueBuildID (CCTK_POINTER_TO_CONST IN cctkGH) +USES FUNCTION UniqueBuildID + +# Return a pointer to an unmodifiable C string +# which contains a unique ID for this simulation +CCTK_POINTER_TO_CONST \ +FUNCTION UniqueSimulationID (CCTK_POINTER_TO_CONST IN cctkGH) +USES FUNCTION UniqueSimulationID + +# Return a pointer to an unmodifiable C string +# which contains a unique ID for this run +CCTK_POINTER_TO_CONST \ +FUNCTION UniqueRunID (CCTK_POINTER_TO_CONST IN cctkGH) +USES FUNCTION UniqueRunID + + + +# Check whether existing output files should be truncated or not +CCTK_INT FUNCTION IO_TruncateOutputFiles (CCTK_POINTER_TO_CONST IN cctkGH) +REQUIRES FUNCTION IO_TruncateOutputFiles + + + +# The setup of the system +CCTK_INT FUNCTION MultiPatch_GetSystemSpecification \ + (CCTK_INT OUT maps) +USES FUNCTION MultiPatch_GetSystemSpecification diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch.par new file mode 100644 index 000000000..71550e788 --- /dev/null +++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch.par @@ -0,0 +1,45 @@ +ActiveThorns = " + Boundary + Carpet + CarpetIOF5 + CarpetLib + CarpetRegrid2 + CartGrid3D + CoordBase + Coordinates + 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 = 0 + +IO::out_dir = $parfile + +Carpet::domain_from_multipatch = yes +CartGrid3D::type = "multipatch" +Coordinates::coordinate_system = "Thornburg04" +Coordinates::h_cartesian = 0.25 +Coordinates::h_radial = 0.25 +Coordinates::sphere_inner_radius = 1.0 +Coordinates::sphere_outer_radius = 3.0 +Coordinates::n_angular = 8 + +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" diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined.par b/CarpetDev/CarpetIOF5/par/iof5-refined.par new file mode 100644 index 000000000..aa27521b4 --- /dev/null +++ b/CarpetDev/CarpetIOF5/par/iof5-refined.par @@ -0,0 +1,48 @@ +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 = 0 + +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" diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform.par b/CarpetDev/CarpetIOF5/par/iof5-uniform.par new file mode 100644 index 000000000..87e758370 --- /dev/null +++ b/CarpetDev/CarpetIOF5/par/iof5-uniform.par @@ -0,0 +1,39 @@ +ActiveThorns = " + Boundary + Carpet + CarpetIOF5 + CarpetLib + 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 = 2 + +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 diff --git a/CarpetDev/CarpetIOF5/param.ccl b/CarpetDev/CarpetIOF5/param.ccl new file mode 100644 index 000000000..f9da760a4 --- /dev/null +++ b/CarpetDev/CarpetIOF5/param.ccl @@ -0,0 +1,66 @@ +# Parameter definitions for thorn CarpetIOF5 + +SHARES: IO + +USES STRING out_dir AS IO_out_dir + +USES INT out_timesteps_per_file + + + +PRIVATE: + +STRING out_dir "Output directory (overrides IO::out_dir)" STEERABLE=always +{ + "^$" :: "Empty: use IO::out_dir" + ".+" :: "Not empty: directory name" +} "" + +KEYWORD file_content "Create one file for every x" STEERABLE=always +{ + "group" :: "" + "thorn" :: "" + "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 +{ + "" :: "" +} "output" + +INT processor_digits "Minimum number of digits for processor number" STEERABLE=always +{ + 0:* :: "" +} 6 + +STRING out_extension "File name extension" STEERABLE=always +{ + "" :: "File extension (including a leading dot, if desired)" +} ".f5" + +BOOLEAN create_subdirs "Create subdirectories for the output files to reduce the number of files per directory" STEERABLE=always +{ +} "no" + +BOOLEAN one_dir_per_file "Create one subdirectory per output file to reduce locking overhead" STEERABLE=always +{ +} "no" + + + +BOOLEAN output_symmetry_points "Output symmetry and inter-patch boundary points" STEERABLE=always +{ +} "no" + +BOOLEAN output_ghost_points "Output ghost points" STEERABLE=always +{ +} "no" + +BOOLEAN output_boundary_points "Output outer boundary points" STEERABLE=always +{ +} "yes" diff --git a/CarpetDev/CarpetIOF5/schedule.ccl b/CarpetDev/CarpetIOF5/schedule.ccl new file mode 100644 index 000000000..494ce02fe --- /dev/null +++ b/CarpetDev/CarpetIOF5/schedule.ccl @@ -0,0 +1,27 @@ +# Schedule definitions for thorn CarpetIOF5 + +SCHEDULE F5_Output AT initial +{ + LANG: C + OPTIONS: global-late +} "Create an output file" + + + +SCHEDULE F5_Poison AT initial AFTER F5_Output BEFORE F5_Input +{ + LANG: C + OPTIONS: global-late +} "Poison all variables" + +SCHEDULE F5_Input AT initial AFTER F5_Output +{ + LANG: C + OPTIONS: global-late +} "Read from file" + +SCHEDULE F5_Check AT initial AFTER F5_Input +{ + LANG: C + OPTIONS: global-late +} "Check all variables for poison" diff --git a/CarpetDev/CarpetIOF5/src/attributes.cc b/CarpetDev/CarpetIOF5/src/attributes.cc new file mode 100644 index 000000000..9ca1f109c --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/attributes.cc @@ -0,0 +1,184 @@ +#include + +#include +#include +#include +#include + +#include + +#include + +#include "iof5.hh" + + + +namespace CarpetIOF5 { + + using namespace std; + + + + // Write an int attribute + 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 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)); + + return error_flag; + } + + + + // Write a double attribute + 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 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)); + + return error_flag; + } + + + + // Write a string attribute + 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 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)); + + return error_flag; + } + + + + // Write an array of int attributes + 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 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)); + + return error_flag; + } + + + + // Write an array of double attributes + 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 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)); + + return error_flag; + } + + + + // 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 error_flag = false; + + size_t maxstrlen = 0; + for (int i=0; i svalue (nvalues * (maxstrlen+1)); + for (int i=0; i +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include "iof5.hh" + + + +namespace CarpetIOF5 { + + using namespace std; + using namespace Carpet; + + + + // Use a class for reading data, so that we have an easy way to pass + // variables between the various iterators + class input_iterator_t { + cGH *cctkGH; + + double time; + char const *gridname; + char const *topologyname; + int index_depth; // 0=vertex, 1=cell + int topological_dimension; + char const *fieldname; + char const *fragmentname; + + // string chartname; + + public: + input_iterator_t (cGH *const cctkGH_) + : cctkGH(cctkGH_), + time(nan), + gridname(NULL), + topologyname(NULL), index_depth(-1), topological_dimension(-1), + fieldname(NULL), + fragmentname(NULL) + { + } + + + + private: + + 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=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; + 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; + } + + 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); + } + + void read_topology (F5Path *const path) + { + indent_t indent; + + herr_t herr; + + cout << indent << "topology=\"" << topologyname << "\"" + << " (" << (index_depth==0 ? "vertex" : "cell") << ")\n"; + + // 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); + if (stat.type == H5G_LINK) { + char linkval[100000]; + herr = H5Gget_linkval + (path->Grid_hid, topologyname, sizeof linkval, linkval); + assert (not herr); + indent_t indent; + cout << indent << "alias for topology \"" << linkval << "\"\n" + << indent << "ignoring this topology\n"; + return; + } + + // F5iterate_topology_fields + // (path, NULL, field_iterator, this, chartname.c_str(), NULL); + F5iterate_topology_fields (path, NULL, field_iterator, this, NULL, NULL); + } + + void read_field (F5Path *const path) + { + indent_t indent; + cout << indent << "field=\"" << fieldname << "\"\n"; + + if (strcmp(fieldname, "GRID::COORDINATES")!=0 and + strcmp(fieldname, "GRID::r")!=0) + { + indent_t indent; + cout << indent << "ignoring this field\n"; + return; + } + + F5iterate_field_fragments (path, NULL, fragment_iterator, this); + } + + void read_fragment (F5Path *const path) + { + indent_t indent; + + int iret; + void *pret; + + cout << indent << "fragment=\"" << fragmentname << "\"\n"; + + + + // Determine refinement level from the topology + // (Could also use topology name instead) + hsize_t hreffact[FIBER_MAX_RANK]; + iret = F5LAget_dimensions(path->Topology_hid, + FIBER_HDF5_REFINEMENT_INFO, hreffact); + assert (iret == dim); + hsize_t hreffact2[FIBER_MAX_RANK]; + pret = F5Tpermute_dimensions(path, dim, hreffact2, hreffact); + assert (pret); + ivect const reffact = h2v(hreffact2); + int rl; + for (rl=0; rl= 0); + + cout << indent << "dataset=\"" << name << "\"\n"; + + assert (var >= 0); + { + char *const fullname = CCTK_FullName(var); + cout << indent << "variable=" << fullname << "\n"; + free (fullname); + } + + + + // Determine fragment properties + H5O_info_t info; + herr = + H5Oget_info_by_name (path->Field_hid, fragmentname, &info, H5P_DEFAULT); + assert (not herr); + bool const fragment_is_group = info.type == H5O_TYPE_GROUP; + hid_t field; + if (fragment_is_group) { + field = H5Gopen (path->Field_hid, fragmentname, H5P_DEFAULT); + assert (field >= 0); + } else { + field = path->Field_hid; + } + + hid_t const fragment = H5Dopen(field, name, H5P_DEFAULT); + assert (fragment); + + // Check index depth + int index_depth_; + iret = F5Tget_index_depth(path, &index_depth_); + assert (iret); + assert (index_depth_ == index_depth); + + // Read the fragment offset. This is stored with the dataset + // group. + hsize_t hoff[FIBER_MAX_RANK]; + iret = F5LAget_dimensions(fragment_is_group ? field : fragment, + FIBER_FRAGMENT_OFFSET_ATTRIBUTE, hoff); + assert (iret == dim); + hsize_t hoff2[FIBER_MAX_RANK]; + pret = F5Tpermute_dimensions(path, dim, hoff2, hoff); + assert (pret); + ivect const foff = h2v(hoff2); + assert (all(foff>=0)); + +#if 0 + // Read the fragment size. This is stored with the field -- why + // is this different from the offset? + hsize_t hlen[FIBER_MAX_RANK]; + iret = + F5LAget_dimensions(path->Field_hid, + FIBER_FIELD_DATASPACE_DIMENSIONS_ATTRIBUTE, hlen); + assert (iret == dim); + hsize_t hlen2[FIBER_MAX_RANK]; + pret = F5Tpermute_dimensions(path, dim, hlen2, hlen); + assert (pret); + ivect const flen = h2v(hlen2); + assert (all(flen>=0)); +#endif + hid_t const space = H5Dget_space(fragment); + assert (space>=0); + iret = H5Sget_simple_extent_ndims(space); + 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); + ivect const flen = h2v(hlen2); + assert (all(flen>=0)); + herr = H5Sclose (space); + assert (not herr); + + ibbox const fbox (foff, foff+flen-1, 1); + { + indent_t indent; + cout << indent << "dataset bbox is " << foff << ":" << foff+flen << "\n"; + } + + + + // 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 is too expensive; only traverse Carpet's data + // structures instead. (This should be implemented via a + // gh::locate_positions, which will turn require support from + // the internal tree data structure.) + BEGIN_COMPONENT_LOOP(cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + + ivect lbnd, lsh, lghosts, ughosts; + ivect imin, imax, ioff, ilen; + for (int d=0; dprocessor(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 (fragment, H5T_NATIVE_DOUBLE, memspace, filespace, + H5P_DEFAULT, data); + assert (not herr); + + herr = H5Sclose(memspace); + assert (not herr); + herr = H5Sclose(filespace); + assert (not herr); + + } END_COMPONENT_LOOP; + + + + herr = H5Dclose(fragment); + assert (not herr); + + if (fragment_is_group) { + herr = H5Gclose(field); + assert (not herr); + } + } + + + + public: + + void iterate (hid_t const object) + { + F5iterate_timeslices (object, NULL, timeslice_iterator, this); + } + + static + 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); + return 0; + } + + static + 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); + 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) + { + 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); + return 0; + } + + static + 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); + return 0; + } + + static + 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); + return 0; + } + + }; // class input_iterator_t + + + + extern "C" + void F5_Input (CCTK_ARGUMENTS) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + herr_t herr; + + + + assert (is_global_mode()); + CCTK_VInfo (CCTK_THORNSTRING, "F5_Input: iteration=%d", cctk_iteration); + + + + // Open file + string const basename = + generate_basename (cctkGH, CCTK_VarIndex("grid::r")); + int const myproc = CCTK_MyProc(cctkGH); + int const nprocs = CCTK_nProcs(cctkGH); + for (int proc=myproc; ; proc+=nprocs) { + string const name = + create_filename (cctkGH, basename, proc, false); + + bool file_exists; + H5E_BEGIN_TRY { + file_exists = H5Fis_hdf5(name.c_str()) > 0; + } H5E_END_TRY; + if (not file_exists) break; + + indent_t indent; + cout << indent << "process=" << proc << "\n"; + + hid_t const file = H5Fopen (name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + assert (file >= 0); + + // Iterate over all time slices + input_iterator_t iterator(cctkGH); + iterator.iterate (file); + + // Close file + herr = H5Fclose (file); + assert (not herr); + } + } + +} // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc.old b/CarpetDev/CarpetIOF5/src/iof5.cc.old new file mode 100644 index 000000000..fe780030b --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/iof5.cc.old @@ -0,0 +1,643 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + + + +namespace CarpetIOF5 { + + using namespace std; + using namespace Carpet; + + + + // File mode for creating directories + int const mode = 0755; + + // A nan + CCTK_REAL const nan = numeric_limits::quiet_NaN(); + + // Indentation + inline string indent (int const n) + { + int const width = 3; + return string(width*n, ' '); + } + + + + typedef vect hvect; + typedef vect fvect; + + + + static inline + hvect v2h (ivect const& a) + { + return hvect(a); + } + + static inline + F5_vec3_point_t v2p (rvect const& a) + { + F5_vec3_point_t res; + res.x = a[0]; + res.y = a[1]; + res.z = a[2]; + return res; + } + + static inline + F5_vec3_float_t v2f (rvect const& a) + { + F5_vec3_float_t res; + res.x = a[0]; + res.y = a[1]; + res.z = a[2]; + return res; + } + + static inline + F5_vec3_double_t v2d (rvect const& a) + { + F5_vec3_double_t res; + res.x = a[0]; + res.y = a[1]; + res.z = a[2]; + return res; + } + + + + // Generate a good grid name (simulation name) + string + generate_gridname (cGH const *const cctkGH) + { +#if 0 + char const *const gridname = (char const*) UniqueSimulationID(cctkGH); + 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); + char *const p = strstr(parfilename, ".par"); + if (p) *p = '\0'; + char *const s = strrchr(parfilename, '/'); + char *const gridname = s ? s+1 : parfilename; + return gridname; + } + + + + // Generate a good file name ("alias") for a variable + string + generate_basename (cGH const *const cctkGH, + int const variable) + { + DECLARE_CCTK_PARAMETERS; + + assert (variable >= 0); + + ostringstream filename_buf; + + if (CCTK_EQUALS (file_content, "variable")) { + char *const varname = CCTK_FullName(variable); + assert (varname); + for (char *p = varname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << varname; + free (varname); + } + else if (CCTK_EQUALS (file_content, "group")) { + char *const groupname = CCTK_GroupNameFromVarI(variable); + assert (groupname); + for (char *p = groupname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << groupname; + free (groupname); + } + else if (CCTK_EQUALS (file_content, "thorn")) { + char const *const impname = CCTK_ImpFromVarI(variable); + char *const thornname = strdup(impname); + assert (thornname); + char *const colon = strchr(thornname, ':'); + assert (colon); + *colon = '\0'; + for (char *p = thornname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << thornname; + free (thornname); + } + else if (CCTK_EQUALS (file_content, "everything")) { + filename_buf << out_filename; + } + else { + 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; + } + + return filename_buf.str(); + } + + + + // Create the final file name on a particular processor + string + create_filename (cGH const *const cctkGH, + string const basename, + bool const create_directories) + { + DECLARE_CCTK_PARAMETERS; + + int const proc = CCTK_MyProc(cctkGH); + int const nprocs = CCTK_nProcs(cctkGH); + + bool const use_IO_out_dir = strcmp(out_dir, "") == 0; + string path = use_IO_out_dir ? IO_out_dir : out_dir; + + if (create_subdirs) { + if (nprocs >= 10000) { + ostringstream buf; + buf << path + "/" + << basename + << ".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); + } + } + if (nprocs >= 100) { + ostringstream buf; + buf << path + "/" + << basename + << ".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); + } + } + } + if (one_dir_per_file and nprocs >= 1) { + ostringstream buf; + buf << path + << basename + << ".p" << setw (processor_digits) << setfill ('0') + << proc + << "/"; + path = buf.str(); + if (create_directories) { + 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(); + } + + + + extern "C" + void F5_Output (CCTK_ARGUMENTS) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + herr_t herr; + + + + assert (is_global_mode()); + CCTK_VInfo (CCTK_THORNSTRING, "F5_Output: iteration=%d", cctk_iteration); + + string const gridname = generate_gridname(cctkGH); + + + + // Open file + static bool first_time = true; + + string const basename = + generate_basename (cctkGH, CCTK_VarIndex("grid::r")); + string const name = + create_filename (cctkGH, basename, first_time); + + 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); + first_time = false; + + + + BEGIN_REFLEVEL_LOOP (cctkGH) { + DECLARE_CCTK_ARGUMENTS; + +#if 0 + // Write data +#warning "TODO: Save the coordinates in double precision" + F5Path *const path = + F5Fwrite_uniform_cartesian3D + (file, cctk_time, gridname, &v2p(lower), &v2f(delta), &v2h(lsh)[0], + "radius" /* group name */, H5T_NATIVE_DOUBLE, r, NULL, F5P_DEFAULT); +#endif + + // Define grid hierarchy + ivect const reffact = spacereffacts.AT(reflevel); + F5Path *const path = F5Rcreate_vertex_refinement3D + (file, cctk_time, gridname.c_str(), &v2h(reffact)[0], NULL); + + // Define iteration + F5Rset_timestep (path, cctk_iteration); + +#if 0 + F5Path *refpath = NULL; + static CCTK_REAL last_root_time = -1; + if (reflevel == 0) { + last_root_time = cctk_time; + } else { + refpath = F5Rcreate_relative_vertex_Qrefinement3D + (file, cctk_time, gridname, &v2h(reffact)[0], + last_root_time, &v2h(ivect(1))[0]); + } +#endif + + + + BEGIN_LOCAL_MAP_LOOP (cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + + // Determine level coordinates + ivect gsh; + rvect origin, delta; + rvect lower, upper; + for (int d=0; dGrid_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); + cout << indent(4) << "alias for topology \"" << linkval << "\"\n" + << indent(4) << "ignoring this topology\n"; + return; + } + + F5iterate_topology_fields (path, NULL, field_iterator, this, NULL, NULL); + } + + void read_field (F5Path *const path) + { + cout << indent(4) << "field=\"" << fieldname << "\"\n"; + F5iterate_field_fragments (path, NULL, fragment_iterator, this); + } + + void read_fragment (F5Path *const path) + { + herr_t herr; + + cout << indent(5) << "fragment=\"" << fragmentname << "\"\n"; + + if (strcmp(fieldname, "radius") == 0) { + + hid_t const group = path->Field_hid; + read_variable (group, fragmentname, CCTK_VarIndex("grid::r")); + + } else if (strcmp(fieldname, "coordinates") == 0) { + + hid_t const group = + H5Gopen (path->Field_hid, fragmentname, H5P_DEFAULT); + assert (group); + read_variable (group, "x", CCTK_VarIndex("grid::x")); + read_variable (group, "y", CCTK_VarIndex("grid::y")); + read_variable (group, "z", CCTK_VarIndex("grid::z")); + herr = H5Gclose (group); + assert (not herr); + + } else { + // do nothing + } + } + + void read_variable (hid_t const group, char const *const name, + int const var) + { + assert (group >= 0); + assert (name); + assert (var >= 0); + + cout << indent(6) << "dataset=\"" << name << "\"\n"; + // hid_t const fragment = H5Dopen(group, name); + // assert (fragment); + // + // hid_t const space = H5Dget_space(fragment); + // assert (space); + // int const rank = H5Sget_simple_extent_dims(space, NULL, NULL); + // assert (rank>=0); + // + // assert (rank == cctk_dim); + // + // H5Sclose(space_id); + } + + + + public: + + void iterate (hid_t const object) + { + F5iterate_timeslices (object, NULL, timeslice_iterator, this); + } + + static + herr_t timeslice_iterator (F5Path *const path, double const time, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->time = time; + iterator->read_timeslice (path); + return 0; + } + + static + herr_t grid_iterator (F5Path *const path, char const *const gridname, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->gridname = gridname; + 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) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->topologyname = topologyname; + iterator->index_depth = index_depth; + iterator->topological_dimension = topological_dimension; + iterator->read_topology (path); + return 0; + } + + static + herr_t field_iterator (F5Path *const path, char const *const fieldname, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->fieldname = fieldname; + iterator->read_field (path); + return 0; + } + + static + herr_t fragment_iterator (F5Path *const path, + char const *const fragmentname, + void *const userdata) + { + iterator_t* const iterator = (iterator_t*)userdata; + iterator->fragmentname = fragmentname; + iterator->read_fragment (path); + return 0; + } + + }; + + + + extern "C" + void F5_Input (CCTK_ARGUMENTS) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + herr_t herr; + + + + assert (is_global_mode()); + CCTK_VInfo (CCTK_THORNSTRING, "F5_Input: iteration=%d", cctk_iteration); + + + + // Open file + string const basename = + generate_basename (cctkGH, CCTK_VarIndex("grid::r")); + string const name = + create_filename (cctkGH, basename, false); + + hid_t const file = H5Fopen (name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); + assert (file >= 0); + + // Iterate over all time slices + iterator_t iterator(cctkGH); + iterator.iterate (file); + + // Close file + herr = H5Fclose (file); + assert (not herr); + } + +} // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/iof5.hh b/CarpetDev/CarpetIOF5/src/iof5.hh new file mode 100644 index 000000000..71524fa2c --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/iof5.hh @@ -0,0 +1,214 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + + + +// This requires defining a boolean local variable "error_flag", +// initialised to false +#define FAILWARN(_expr) \ + failwarn(error_flag, _expr, __LINE__, __FILE__, CCTK_THORNSTRING, #_expr) + +template +static +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); + error_flag = true; + } + return expr; +} + + + +namespace CarpetIOF5 { + class indent_t; +}; +ostream& operator<<(ostream& os, CarpetIOF5::indent_t const& indent); + + + +namespace CarpetIOF5 { + + using namespace std; + using namespace Carpet; + + + + // Indentation + class indent_t { + static int const width = 3; + static int level; + public: + indent_t() { ++level; } + ~indent_t() { --level; } + ostream& output(ostream& os) const; + }; + + + + // File mode for creating directories + int const mode = 0755; + + // A nan + CCTK_REAL const nan = numeric_limits::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; + + + + // Data types for HDF5 and F5 types + typedef vect hvect; + + // Conversion operators for these datatypes + static inline + hvect v2h (ivect const& a) + { + return hvect(a); + } + static inline + ivect h2v (hvect const& a) + { + return ivect(a); + } + static inline + ivect h2v (hsize_t const* const a) + { + ivect res; + for (int d=0; d +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h" + +#include +#include + +#include + +#include "iof5.hh" + + + +namespace CarpetIOF5 { + + using namespace std; + using namespace Carpet; + + + + class output_iterator_t { + cGH *const cctkGH; + + bool const is_multipatch; + + string gridname; + string chartname; + string fragmentname; + string topologyname; + + public: + output_iterator_t (cGH *const cctkGH_) + : cctkGH(cctkGH_), + is_multipatch + (CCTK_IsFunctionAliased("MultiPatch_GetSystemSpecification")) + { + } + + void iterate (hid_t const file) + { + int const rl = 0; + ENTER_LEVEL_MODE(cctkGH, rl) { + output_simulation (file); + } LEAVE_LEVEL_MODE; + } + + + + private: + + void output_simulation (hid_t const file) + { + assert (is_level_mode() and reflevel==0); + DECLARE_CCTK_ARGUMENTS; + indent_t indent; + + gridname = generate_gridname(cctkGH); + chartname = generate_chartname(cctkGH); + + cout << indent << "simulation=" << gridname << "\n"; + + ivect const reffact = spacereffacts.AT(reflevel); + F5Path *const globalpath = F5Rcreate_vertex_refinement3D + (file, cctk_time, gridname.c_str(), &v2h(reffact)[0], + chartname.c_str()); + + // Choose (arbitrarily) the root level as default topology, for + // readers which don't understand AMR + F5Rlink_default_vertex_topology (globalpath, &v2h(reffact)[0]); + + // Define iteration + F5Rset_timestep (globalpath, cctk_iteration); + + // Attach Cactus/Carpet metadata + write_metadata (cctkGH, globalpath->Grid_hid); + + // Close topology + F5close (globalpath); + + // Iterate over all maps (on the root level) + BEGIN_MAP_LOOP(cctkGH, CCTK_GF) { + output_map (file); + } END_MAP_LOOP; + } + + + + void output_map (hid_t const file) + { + DECLARE_CCTK_ARGUMENTS; + indent_t indent; + + herr_t herr; + + assert (is_singlemap_mode() and reflevel==0); + + cout << indent << "map=" << Carpet::map << "\n"; + + fragmentname = generate_fragmentname(cctkGH, Carpet::map); + + // Iterate over all refinement levels of this map + int const m = Carpet::map; + BEGIN_GLOBAL_MODE(cctkGH) { + BEGIN_REFLEVEL_LOOP(cctkGH) { + if (vhh.AT(m)->local_components(reflevel) > 0) { + // Continue only if this process owns a component of this + // map (on this level) + ENTER_SINGLEMAP_MODE(cctkGH, m, CCTK_GF) { + output_reflevel (file); + } LEAVE_SINGLEMAP_MODE; + } + } END_REFLEVEL_LOOP; + } END_GLOBAL_MODE; + } + + + + void output_reflevel (hid_t const file) + { + DECLARE_CCTK_ARGUMENTS; + indent_t indent; + + cout << indent << "reflevel=" << reflevel << "\n"; + + ivect const reffact = spacereffacts.AT(reflevel); + topologyname = generate_topologyname(cctkGH, Carpet::map, reffact); + + // Define grid hierarchy + int const indexdepth = 0; // vertices + F5Path *const path = + F5Rcreate_coordinate_topology (file, &cctk_time, + gridname.c_str(), chartname.c_str(), + topologyname.c_str(), + indexdepth, + dim, dim, &v2h(reffact)[0]); + + // Determine level coordinates + ivect gsh; + rvect origin, delta; + rvect lower, upper; + for (int d=0; d= 0); + { + char *const fullname = CCTK_FullName(var); + cout << indent << "variable=" << fullname << "\n"; + free (fullname); + } + + int const group = CCTK_GroupIndexFromVarI(var); + int const v0 = CCTK_FirstVarIndexI(group); + + 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 == cctk_dim); + + int const timelevel = 0; + + // Determine component coordinates + int const dim = cctk_dim; + ivect gsh; + ivect lbnd, lsh, lghosts, ughosts; + ivect imin, imax, ioff, ilen; + ivect const izero(0); + for (int d=0; d= 0); + + if (group == coordinates_group) { + + // Special case + if (var >= v0 and var < v0+dim) { + tensortype = tt_vector; + } else if (var == v0+dim) { + tensortype = tt_scalar; + } else { + assert (0); + } + + } else { + + // TODO: use the tensor type instead + switch (groupdata.numvars) { + case 1: + tensortype = tt_scalar; break; + case 3: + tensortype = tt_vector; break; + default: + // fallback: output group as scalars + tensortype = tt_scalar; break; + } + + } + + if (tensortype != tt_scalar and var != v0) { + // TODO: don't do this; instead, keep track of which groups + // have been output + char *const fullname = CCTK_FullName(var); + indent_t indent; + cout << indent << "skipping output since it is not the first variable in its group\n"; + free (fullname); + return; + } + + + + switch (tensortype) { + + case tt_scalar: { + // Scalar field + CCTK_REAL const *const rdata = + (CCTK_REAL const*)CCTK_VarDataPtrI (cctkGH, timelevel, var); + assert (rdata); + + int const vartype = CCTK_VarTypeI(var); + assert (vartype == CCTK_VARIABLE_REAL); + vector idata (prod(ilen)); + for (int k=0; k idata[cctk_dim]; + void const* data[cctk_dim]; + for (int d=0; d= 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 + + // Create a group to hold all metadata + hid_t const group = + H5Gcreate (file, metadata_group, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert (group >= 0); + + // General information + WriteAttribute (group, "Cactus version", CCTK_FullVersion()); + + // Unique identifiers + if (CCTK_IsFunctionAliased ("UniqueConfigID")) { + WriteAttribute + (group, "config id", (char const*) UniqueConfigID(cctkGH)); + } + if (CCTK_IsFunctionAliased ("UniqueBuildID")) { + WriteAttribute + (group, "build id", (char const*) UniqueBuildID(cctkGH)); + } + if (CCTK_IsFunctionAliased ("UniqueSimulationID")) { + WriteAttribute + (group, "simulation id", (char const*) UniqueSimulationID(cctkGH)); + } + if (CCTK_IsFunctionAliased ("UniqueRunID")) { + WriteAttribute + (group, "run id", (char const*) UniqueRunID(cctkGH)); + } + + // 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); + + // 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); + + // Grid structure + string const gs = serialise_grid_structure (cctkGH); + WriteLargeAttribute (group, grid_structure, gs.c_str()); + + herr = H5Gclose (group); + assert (not herr); + } + + + + extern "C" + void F5_Output (CCTK_ARGUMENTS) + { + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + herr_t herr; + + + + assert (is_global_mode()); + CCTK_VInfo (CCTK_THORNSTRING, "F5_Output: iteration=%d", cctk_iteration); + + + + // We don't know how to open multiple files yet + assert (strcmp (file_content, "everything") == 0); + + // Open file + static bool first_time = true; + + string const basename = + generate_basename (cctkGH, CCTK_VarIndex("grid::r")); + int const myproc = CCTK_MyProc(cctkGH); + int const proc = myproc; + string const name = + create_filename (cctkGH, basename, proc, first_time); + + indent_t indent; + cout << indent << "process=" << proc << "\n"; + + 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); + first_time = false; + + // write_metadata (cctkGH, file); + + output_iterator_t iterator(cctkGH); + iterator.iterate (file); + + // Close file + herr = H5Fclose (file); + assert (not herr); + } + +} // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/poison.cc b/CarpetDev/CarpetIOF5/src/poison.cc new file mode 100644 index 000000000..6029b55a8 --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/poison.cc @@ -0,0 +1,77 @@ +#include +#include +#include + +#include +#include + +#include "iof5.hh" + + + +namespace CarpetIOF5 { + + using namespace std; + using namespace Carpet; + + + + extern "C" + void F5_Poison (CCTK_ARGUMENTS) + { + DECLARE_CCTK_PARAMETERS; + + 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) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + LC_LOOP3(F5_Poison, i,j,k, + 0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k); + x[ind3d] = nan; + y[ind3d] = nan; + z[ind3d] = nan; + r[ind3d] = nan; + } LC_ENDLOOP3(F5_Poison); + } END_LOCAL_COMPONENT_LOOP; + } END_LOCAL_MAP_LOOP; + } END_REFLEVEL_LOOP; + } + + + + extern "C" + void F5_Check (CCTK_ARGUMENTS) + { + DECLARE_CCTK_PARAMETERS; + + 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) { + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + LC_LOOP3(F5_Check, i,j,k, + 0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + 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])); + } LC_ENDLOOP3(F5_Check); + } END_LOCAL_COMPONENT_LOOP; + } END_LOCAL_MAP_LOOP; + } END_REFLEVEL_LOOP; + } + +} // end namespace CarpetIOF5 diff --git a/CarpetDev/CarpetIOF5/src/util.cc b/CarpetDev/CarpetIOF5/src/util.cc new file mode 100644 index 000000000..57ccfd6db --- /dev/null +++ b/CarpetDev/CarpetIOF5/src/util.cc @@ -0,0 +1,282 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include "iof5.hh" + + + +ostream& operator<<(ostream& os, CarpetIOF5::indent_t const& indent) +{ return indent.output(os); } + + + +namespace CarpetIOF5 { + + using namespace std; + using namespace Carpet; + + + + // Indentation + int indent_t::level = 0; + + ostream& indent_t::output(ostream& os) const + { + return os << string(width*level, ' '); + } + + + + // Generate a good file name ("alias") for a variable + string + generate_basename (cGH const *const cctkGH, + int const variable) + { + DECLARE_CCTK_PARAMETERS; + + ostringstream filename_buf; + + if (CCTK_EQUALS (file_content, "variable")) { + assert (variable >= 0); + char *const varname = CCTK_FullName(variable); + assert (varname); + for (char *p = varname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << varname; + free (varname); + } + else if (CCTK_EQUALS (file_content, "group")) { + assert (variable >= 0); + char *const groupname = CCTK_GroupNameFromVarI(variable); + assert (groupname); + for (char *p = groupname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << groupname; + free (groupname); + } + else if (CCTK_EQUALS (file_content, "thorn")) { + assert (variable >= 0); + char const *const impname = CCTK_ImpFromVarI(variable); + char *const thornname = strdup(impname); + assert (thornname); + char *const colon = strchr(thornname, ':'); + assert (colon); + *colon = '\0'; + for (char *p = thornname; *p; ++p) { + *p = tolower(*p); + } + filename_buf << thornname; + free (thornname); + } + else if (CCTK_EQUALS (file_content, "everything")) { + filename_buf << out_filename; + } + else { + 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; + } + + return filename_buf.str(); + } + + + + // Create the final file name on a particular processor + string + create_filename (cGH const *const cctkGH, + string const basename, + int const proc, + bool const create_directories) + { + DECLARE_CCTK_PARAMETERS; + + bool const use_IO_out_dir = strcmp(out_dir, "") == 0; + string path = use_IO_out_dir ? IO_out_dir : out_dir; + + if (create_subdirs) { + { + ostringstream buf; + buf << path + "/" + << basename + << ".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); + } + } + { + ostringstream buf; + buf << path + "/" + << basename + << ".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); + } + } + } + if (one_dir_per_file) { + ostringstream buf; + buf << path + << basename + << ".p" << setw (processor_digits) << setfill ('0') + << proc + << "/"; + path = buf.str(); + if (create_directories) { + 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(); + } + + + + // Generate a good grid name (simulation name) + string + generate_gridname (cGH const *const cctkGH) + { +#if 0 + char const *const gridname = (char const*) UniqueSimulationID(cctkGH); + 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); + char const *const suffix = ".par"; + if (strlen(parfilename) >= strlen(suffix)) { + char *const p = parfilename + strlen(parfilename) - strlen(suffix); + if (strcmp(p, suffix) == 0) { + *p = '\0'; + } + } + char *const s = strrchr(parfilename, '/'); + char *const p = s ? s+1 : parfilename; + + return p; + } + + // Generate a good tensor basis name (chart name) + string + generate_chartname (cGH const *const cctkGH) + { + // Assume a global tensor basis, where all maps use the same chart + return FIBER_HDF5_DEFAULT_CHART; + +#if 0 + // TODO: Add some interesting information here, e.g. the name of + // the multi-patch system, or a good name of the patch? + ostringstream buf; + buf << "Map_" << m; + return buf.str(); +#endif + } + + // Generate a good fragment name (map name) + string + generate_fragmentname (cGH const *const cctkGH, int const m) + { + ostringstream buf; + buf << "Map_" << m; + return buf.str(); + } + + int + interpret_fragmentname (cGH const *const cctkGH, + char const *const fragmentname) + { + int m = -1; + sscanf (fragmentname, "Map_%d", &m); + assert (m>=0 and m0) buf << "x"; + buf << reffact[d]; + } + return buf.str(); + } + + + + char const *const grid_structure = "Grid Structure v5"; + + string + serialise_grid_structure (cGH const *const cctkGH) + { + ostringstream buf; + buf << setprecision(17); + buf << grid_structure << "{"; + buf << "maps=" << maps << ","; + buf << "["; + for (int m=0; msuperregions << ","; + // We could omit the regions + buf << "regions:" << vhh.AT(m)->regions.AT(0) << ","; + buf << "ghost_widths:" << vdd.AT(m)->ghost_widths << ","; + buf << "buffer_widths:" << vdd.AT(m)->buffer_widths << ","; + buf << "prolongation_orders_space:" + << vdd.AT(m)->prolongation_orders_space << ","; + buf << "},"; + } + buf << "],"; + buf << "times:" << *tt << ","; + buf << "}."; + return buf.str(); + } + +} // end namespace CarpetIOF5 -- cgit v1.2.3