diff options
Diffstat (limited to 'src/patch/patch_system.cc')
| -rw-r--r-- | src/patch/patch_system.cc | 702 |
1 files changed, 684 insertions, 18 deletions
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc index 98974e6..bbd29f7 100644 --- a/src/patch/patch_system.cc +++ b/src/patch/patch_system.cc @@ -30,6 +30,8 @@ // // patch_system::set_gridfn_to_constant // patch_system::gridfn_copy +// patch_system::ghosted_gridfn_copy +// patch_system::add_to_gridfn // patch_system::add_to_ghosted_gridfn // patch_system::gridfn_norms // patch_system::ghosted_gridfn_norms @@ -39,9 +41,11 @@ // // patch_system::patch_containing_local_xyz // patch_system::radius_in_local_xyz_direction +// patch_system::radii_in_local_xyz_directions // // patch_system::print_unknown_gridfn // patch_system::read_unknown_gridfn +// patch_system::output_unknown_gridfn // // patch_system::synchronize // patch_system::compute_synchronize_Jacobian @@ -51,20 +55,37 @@ /// patch_system::ghost_zone_Jacobian // +#include <assert.h> #include <stdio.h> #include <math.h> #include <string.h> #include <assert.h> #include <limits.h> +#include <sys/stat.h> +#include <sys/types.h> + +#include <vector> +#include <sstream> +#include <string> #include "cctk.h" +#include "CactusBase/IOUtil/src/ioGH.h" #include "config.h" #include "stdc.h" +#include "../gr/gfns.hh" #include "../jtutil/util.hh" #include "../jtutil/array.hh" #include "../jtutil/cpm_map.hh" #include "../jtutil/linear_map.hh" +using jtutil::error_exit; + +#ifdef CCTK_HDF5 +// Some macros to fix compatibility issues as long as 1.8.0 is in beta +// phase +# define H5_USE_16_API 1 +# include <hdf5.h> +#endif #include "coords.hh" #include "grid.hh" @@ -80,7 +101,6 @@ // all the code in this file is inside this namespace namespace AHFinderDirect { -using jtutil::error_exit; //****************************************************************************** //****************************************************************************** @@ -325,12 +345,13 @@ delete[] ghosted_gridfn_storage_; delete[] gridfn_storage_; delete[] ghosted_starting_gpn_; delete[] starting_gpn_; -delete[] all_patches_; for (int pn = N_patches()-1 ; pn >= 0 ; --pn) { delete &ith_patch(pn); } + +delete[] all_patches_; } //****************************************************************************** @@ -342,6 +363,8 @@ delete[] all_patches_; // of this patch system. This function also correctly sets // N_grid_points_ // N_ghosted_grid_points_ +// max_N_additional_points_ +// N_additional_points_ // all_patches_[] // starting_gpn_[] // ghosted_starting_gpn_[] @@ -353,6 +376,9 @@ void patch_system::create_patches(const struct patch_info patch_info_in[], int N_zones_per_right_angle, bool print_msg_flag) { +max_N_additional_points_ = 1; +// N_additional_points_ = 0; +N_additional_points_ = 1; N_grid_points_ = 0; ghosted_N_grid_points_ = 0; for (int pn = 0 ; pn < N_patches() ; ++pn) @@ -461,8 +487,8 @@ const int N_gridfns_in = jtutil::how_many_in_range(min_gfn_in, max_gfn_in); const int ghosted_N_gridfns_in = jtutil::how_many_in_range(ghosted_min_gfn_in, ghosted_max_gfn_in); -const int gfn_stride = N_grid_points(); -const int ghosted_gfn_stride = ghosted_N_grid_points(); +const int gfn_stride = N_grid_points() + max_N_additional_points(); +const int ghosted_gfn_stride = ghosted_N_grid_points() + max_N_additional_points(); const int N_storage = gfn_stride * N_gridfns_in; const int ghosted_N_storage = ghosted_gfn_stride * ghosted_N_gridfns_in; @@ -531,7 +557,7 @@ const patch& plast = ith_patch(N_patches()-1); const fp* const gfn1_first_ptr = & pfirst.gridfn(gfn+1, pfirst.min_irho(), pfirst.min_isigma()); - if (! (gfn1_first_ptr == gfn_last_ptr+1)) + if (! (gfn1_first_ptr == gfn_last_ptr+max_N_additional_points()+1)) then error_exit(PANIC_EXIT, "***** patch_system::setup_gridfn_storage():\n" " nominal-grid gridfns don't partition overall storage array!" @@ -561,7 +587,7 @@ const patch& plast = ith_patch(N_patches()-1); = & pfirst.ghosted_gridfn(ghosted_gfn+1, pfirst.ghosted_min_irho(), pfirst.ghosted_min_isigma()); - if (! (ghosted_gfn1_first_ptr == ghosted_gfn_last_ptr+1)) + if (! (ghosted_gfn1_first_ptr == ghosted_gfn_last_ptr+max_N_additional_points()+1)) then error_exit(PANIC_EXIT, "***** patch_system::setup_gridfn_storage():\n" " ghosted-grid gridfns don't partition overall storage array!" @@ -635,7 +661,7 @@ const patch& plast = ith_patch(N_patches()-1); if (! (p1_first_ptr == p_last_ptr+1)) then error_exit(PANIC_EXIT, "***** patch_system::setup_gridfn_storage():\n" -" nominal-grid patches gridfns don't partition storage for gfn=%d!\n" +" ghosted-grid patches gridfns don't partition storage for gfn=%d!\n" " (this should never happen!)\n" " %s patch (pn=%d) last point at p_last_ptr=%p\n" " %s patch (pn=%d) first point at p1_first_ptr=%p\n" @@ -1489,6 +1515,29 @@ void patch_system::set_gridfn_to_constant(fp a, int dst_gfn) //****************************************************************************** // +// This function sets a ghosted gridfn to a constant value. +// +void patch_system::set_ghosted_gridfn_to_constant(fp a, int ghosted_dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.ghosted_min_irho() ; + irho <= p.ghosted_max_irho() ; ++irho) + { + for (int isigma = p.ghosted_min_isigma() ; + isigma <= p.ghosted_max_isigma() ; + ++isigma) + { + p.ghosted_gridfn(ghosted_dst_gfn, irho,isigma) = a; + } + } + } +} + +//****************************************************************************** + +// // This function copies one (nominal-grid) gridfn to another. // void patch_system::gridfn_copy(int src_gfn, int dst_gfn) @@ -1511,6 +1560,78 @@ void patch_system::gridfn_copy(int src_gfn, int dst_gfn) //****************************************************************************** // +// This function copies one ghosted gridfn to another. +// +void patch_system::ghosted_gridfn_copy(int ghosted_src_gfn, + int ghosted_dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.ghosted_min_irho() ; + irho <= p.ghosted_max_irho() ; ++irho) + { + for (int isigma = p.ghosted_min_isigma() ; + isigma <= p.ghosted_max_isigma() ; + ++isigma) + { + p.ghosted_gridfn(ghosted_dst_gfn, irho,isigma) + = p.ghosted_gridfn(ghosted_src_gfn, irho,isigma); + } + } + } +} + +//****************************************************************************** + +// +// This function copies a ghosted gridfn to a nominal gridfn. +// +void patch_system::ghosted_gridfn_copy_to_nominal(int ghosted_src_gfn, + int dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho) + { + for (int isigma = p.min_isigma() ; isigma <= p.max_isigma() ; + ++isigma) + { + p.gridfn(dst_gfn, irho,isigma) + = p.ghosted_gridfn(ghosted_src_gfn, irho,isigma); + } + } + } +} + +//****************************************************************************** + +// +// This function adds a scalar to a (nominal-grid) gridfn. +// +void patch_system::add_to_gridfn(fp delta, int dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.min_irho() ; + irho <= p.max_irho() ; + ++irho) + { + for (int isigma = p.min_isigma() ; + isigma <= p.max_isigma() ; + ++isigma) + { + p.gridfn(dst_gfn, irho,isigma) += delta; + } + } + } +} + +//****************************************************************************** + +// // This function adds a scalar to a ghosted gridfn. // void patch_system::add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn) @@ -1535,6 +1656,30 @@ void patch_system::add_to_ghosted_gridfn(fp delta, int ghosted_dst_gfn) //****************************************************************************** // +// This function scales a ghosted gridfn with a scalar. +// +void patch_system::scale_ghosted_gridfn(fp factor, int ghosted_dst_gfn) +{ + for (int pn = 0 ; pn < N_patches() ; ++pn) + { + patch& p = ith_patch(pn); + for (int irho = p.ghosted_min_irho() ; + irho <= p.ghosted_max_irho() ; + ++irho) + { + for (int isigma = p.ghosted_min_isigma() ; + isigma <= p.ghosted_max_isigma() ; + ++isigma) + { + p.ghosted_gridfn(ghosted_dst_gfn, irho,isigma) *= factor; + } + } + } +} + +//****************************************************************************** + +// // This function computes norms of a nominal-grid gridfn. // void patch_system::gridfn_norms(int src_gfn, jtutil::norm<fp>& norms) @@ -1757,6 +1902,69 @@ default: return integral; } +fp patch_system::integrate_gridpoint(int unknown_src_gfn, + int ghosted_radius_gfn, + int g_xx_gfn, int g_xy_gfn, int g_xz_gfn, + int g_yy_gfn, int g_yz_gfn, + int g_zz_gfn, + enum patch::integration_method method, + int pn, int irho, int isigma) + const +{ +const patch& p = ith_patch(pn); +const fp val = p.integrate_gridpoint(unknown_src_gfn, + ghosted_radius_gfn, + g_xx_gfn, g_xy_gfn, g_xz_gfn, + g_yy_gfn, g_yz_gfn, + g_zz_gfn, + method, + irho, isigma); + +return val; +} +fp patch_system::integrate_correction(bool src_gfn_is_even_across_xy_plane, + bool src_gfn_is_even_across_xz_plane, + bool src_gfn_is_even_across_yz_plane) + const +{ +fp integral = 1.0; +// +// correct the integral +// for the fact that the patch system may not cover the full 2-sphere +// +switch (type()) + { +case patch_system__full_sphere: + break; +case patch_system__plus_z_hemisphere: + integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0; + break; +case patch_system__plus_xy_quadrant_mirrored: +case patch_system__plus_xy_quadrant_rotating: + integral *= src_gfn_is_even_across_xz_plane ? 2.0 : 0.0; + integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0; + break; +case patch_system__plus_xz_quadrant_mirrored: +case patch_system__plus_xz_quadrant_rotating: + integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0; + integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0; + break; +case patch_system__plus_xyz_octant_mirrored: +case patch_system__plus_xyz_octant_rotating: + integral *= src_gfn_is_even_across_xy_plane ? 2.0 : 0.0; + integral *= src_gfn_is_even_across_xz_plane ? 2.0 : 0.0; + integral *= src_gfn_is_even_across_yz_plane ? 2.0 : 0.0; + break; +default: + error_exit(PANIC_EXIT, +"***** patch_system::integrate_gridfn(): bad patch system type()=(int)%d!\n" +" (this should never happen!)\n", + int(type())); /*NOTREACHED*/ + } + +return integral; +} + //****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -1814,7 +2022,6 @@ else error_exit(ERROR_EXIT, // patch-system symmetries. If the point coincides with the local origin, // we return the dummy value 1.0. // -// Bugs: // Due to the surface-interpolator overhead, repeatedly calling this // function is rather inefficient. // @@ -1940,13 +2147,177 @@ if (status < 0) , status, double(x), double(y), double(z), - p.name(), rho, sigma, + p.name(), double(rho), double(sigma), double(p.drho_of_rho(rho)), double(p.dsigma_of_sigma(sigma))); /*NOTREACHED*/ return xyz_radius; } +struct perpatch { + std::vector<int> n; + std::vector<fp> rho, sigma; + std::vector<fp> radii; + void reserve (int const npoints) + { + n.reserve (npoints); + rho.reserve (npoints); + sigma.reserve (npoints); + radii.reserve (npoints); + } + void addpoint (int const n_, fp const rho_, fp const sigma_) + { + n .push_back (n_); + rho .push_back (rho_); + sigma.push_back (sigma_); + radii.push_back (fp (0.0)); + } +}; + +void patch_system::radii_in_local_xyz_directions(int const ghosted_radius_gfn, + int const npoints, + fp const * const xp, + fp const * const yp, + fp const * const zp, + fp * const radii) + const +{ + std::vector<perpatch> patchpoints (N_patches()); + for (int pn=0; pn<N_patches(); ++pn) { + patchpoints.at(pn).reserve (npoints); + } + + for (int n=0; n<npoints; ++n) { + fp x = xp[n]; + fp y = yp[n]; + fp z = zp[n]; + + if ((x == 0.0) && (y == 0.0) && (z == 0.0)) { + radii[n] = 1.0; + } else { + // + // apply symmetries to map (x,y,z) into that part of the + // 2-sphere which is covered by the patch system + // + switch (type()) + { + case patch_system__full_sphere: + break; + case patch_system__plus_z_hemisphere: + z = fabs(z); + break; + case patch_system__plus_xy_quadrant_mirrored: + case patch_system__plus_xy_quadrant_rotating: + x = fabs(x); + y = fabs(y); + break; + case patch_system__plus_xz_quadrant_mirrored: + case patch_system__plus_xz_quadrant_rotating: + x = fabs(x); + z = fabs(z); + break; + case patch_system__plus_xyz_octant_mirrored: + case patch_system__plus_xyz_octant_rotating: + x = fabs(x); + y = fabs(y); + z = fabs(z); + break; + default: + error_exit(PANIC_EXIT, +"***** patch_system::radii_in_local_xyz_directions():\n" +" unknown patch system type()=(int)%d!\n" +" (this should never happen!)\n", + int(type())); /*NOTREACHED*/ + } + + const patch* p_ptr = patch_containing_local_xyz(x, y, z); + if (p_ptr == NULL) + then error_exit(ERROR_EXIT, +"***** patch_system::radius_in_local_xyz_direction():\n" +" can't find containing patch!\n" +" (this should never happen!)\n" +" [local] (x,y,z)=(%g,%g,%g)\n" + , + double(x), double(y), double(z)); /*NOTREACHED*/ + + const patch& p = *p_ptr; + const fp rho = p. rho_of_xyz(x,y,z); + const fp sigma = p.sigma_of_xyz(x,y,z); + + patchpoints.at(p.patch_number()).addpoint (n, rho, sigma); + } + } // for n + + for (int pn=0; pn<N_patches(); ++pn) { + const patch& p = ith_patch(pn); + perpatch& localpatchpoints = patchpoints.at(pn); + + // Set up the surface interpolator to interpolate the surface + // radius gridfn to the (rho,sigma) coordinates: + // + // Notes on the interpolator setup: + // * The interpolator assumes Fortran subscripting, so we take the + // coordinates in the order (sigma,rho) to match our C + // subscripting in the patch system. + // * To avoid having to set up min/max array subscripts in the + // parameter table, we treat the patch as using 0-origin + // (integer) array subscripts (irho - ghosted_min_irho(), isigma + // - ghosted_min_isigma()). However, we use the usual + // floating-point coordinates. + // + + const int N_dims = 2; + const CCTK_REAL coord_origin[N_dims] + = { p.ghosted_min_sigma(), p.ghosted_min_rho() }; + const CCTK_REAL coord_delta[N_dims] + = { p.delta_sigma(), p.delta_rho() }; + + const int N_interp_points = localpatchpoints.radii.size(); + const int interp_coords_type_code = CCTK_VARIABLE_REAL; + const void* const interp_coords[N_dims] + = { static_cast<const void*>(& localpatchpoints.sigma.front()), + static_cast<const void*>(& localpatchpoints.rho .front()) }; + + const int N_input_arrays = 1; + const CCTK_INT input_array_dims[N_dims] + = { p.ghosted_N_isigma(), p.ghosted_N_irho() }; + const CCTK_INT input_array_type_codes[N_input_arrays] + = { CCTK_VARIABLE_REAL }; + const void* const input_arrays[N_input_arrays] + = { static_cast<const void*> + (p.ghosted_gridfn_data_array(ghosted_radius_gfn)) }; + + const int N_output_arrays = 1; + const CCTK_INT output_array_type_codes[N_output_arrays] + = { CCTK_VARIABLE_REAL }; + void* const output_arrays[N_output_arrays] + = { static_cast<void*>(& localpatchpoints.radii.front()) }; + + const int status = CCTK_InterpLocalUniform + (N_dims, + surface_interp_handle_, + surface_interp_par_table_handle_, + coord_origin, coord_delta, + N_interp_points, + interp_coords_type_code, interp_coords, + N_input_arrays, + input_array_dims, input_array_type_codes, input_arrays, + N_output_arrays, + output_array_type_codes, output_arrays); + if (status < 0) + then error_exit(ERROR_EXIT, +"***** patch_system::radii_in_local_xyz_directions():\n" +" error return (status=%d) from surface interpolator!\n" + , + status); /*NOTREACHED*/ + + for (size_t nn=0; nn<localpatchpoints.radii.size(); ++nn) { + radii[localpatchpoints.n.at(nn)] = localpatchpoints.radii.at(nn); + } + + } // for p +} + //****************************************************************************** //****************************************************************************** //****************************************************************************** @@ -2034,8 +2405,9 @@ fprintf(output_fp, "\n"); const fp dpy = p.dpy_of_rho_sigma(rho, sigma); fprintf(output_fp, "%g\t%g\t%#.15g", - dpx, dpy, p.unknown_gridfn(ghosted_flag, - unknown_gfn, irho,isigma)); + double(dpx), double(dpy), + double(p.unknown_gridfn(ghosted_flag, + unknown_gfn, irho,isigma))); if (print_xyz_flag) then { const fp r = p.unknown_gridfn(radius_is_ghosted_flag, @@ -2049,7 +2421,8 @@ fprintf(output_fp, "\n"); const fp global_z = origin_z() + local_z; fprintf(output_fp, "\t%#.10g\t%#.10g\t%#.10g", - global_x, global_y, global_z); + double(global_x), double(global_y), + double(global_z)); } fprintf(output_fp, "\n"); } @@ -2122,12 +2495,13 @@ int line_number = 1; " dpx=%g dpy=%g\n" , p.name(), unknown_gfn, - irho, p.effective_min_irho(want_ghost_zones), - p.effective_max_irho(want_ghost_zones), - isigma, + int(irho), + p.effective_min_irho(want_ghost_zones), + p.effective_max_irho(want_ghost_zones), + int(isigma), p.effective_min_isigma(want_ghost_zones), p.effective_max_isigma(want_ghost_zones), - dpx, dpy); /*NOTREACHED*/ + double(dpx), double(dpy)); /*NOTREACHED*/ ++line_number; } while ((buffer[0] == '#') || (buffer[0] == '\n')); @@ -2150,8 +2524,8 @@ int line_number = 1; , p.name(), unknown_gfn, line_number, - dpx, dpy, - read_dpx, read_dpy); /*NOTREACHED*/ + double(dpx), double(dpy), + double(read_dpx), double(read_dpy)); /*NOTREACHED*/ p.unknown_gridfn(ghosted_flag, unknown_gfn, irho,isigma) = read_gridfn_value; @@ -2164,6 +2538,298 @@ fclose(input_fp); } //****************************************************************************** + +#ifdef CCTK_HDF5 + +static void WriteAttribute (const hid_t dataset, const char* name, char value); +static void WriteAttribute (const hid_t dataset, const char* name, const char* values); +static void WriteAttribute (const hid_t dataset, const char* name, const char* values, int nvalues); + +static void WriteAttribute (const hid_t dataset, const char* const name, const char value) +{ + WriteAttribute (dataset, name, &value, 1); +} + +static void WriteAttribute (const hid_t dataset, const char* const name, const char* const values) +{ + WriteAttribute (dataset, name, values, strlen(values)); +} + +static void WriteAttribute (const hid_t dataset, const char* const name, const char* const values, const int nvalues) +{ + assert (dataset>=0); + assert (name); + assert (values); + assert (nvalues>=0); + + herr_t herr; + + const hid_t dataspace = H5Screate (H5S_SCALAR); + assert (dataspace>=0); + + const hid_t datatype = H5Tcopy (H5T_C_S1); + assert (datatype>=0); + herr = H5Tset_size (datatype, nvalues); + assert (!herr); + + const hid_t attribute = H5Acreate (dataset, name, datatype, dataspace, H5P_DEFAULT); + assert (attribute>=0); + herr = H5Awrite (attribute, datatype, values); + assert (!herr); + herr = H5Aclose (attribute); + assert (!herr); + + herr = H5Tclose (datatype); + assert (!herr); + + herr = H5Sclose (dataspace); + assert (!herr); +} + +#endif + +// +// This function output an unknown-grid gridfn in HDF5 format to a +// named output file. +// +void patch_system::output_unknown_gridfn + (bool ghosted_flag, int unknown_gfn, + const char gfn_name[], const cGH *cctkGH, + bool output_xyz_flag, bool radius_is_ghosted_flag, + int unknown_radius_gfn, + const char output_file_name[], bool want_ghost_zones) + const +{ +if (want_ghost_zones && !ghosted_flag) + then error_exit(PANIC_EXIT, +"***** patch_system::output_unknown_gridfn(unknown_gfn=%d):\n" +" can't have want_ghost_zones && !ghosted_flag !\n" +, + unknown_gfn); /*NOTREACHED*/ +if (want_ghost_zones && output_xyz_flag && !radius_is_ghosted_flag) + then error_exit(PANIC_EXIT, +"***** patch_system::output_unknown_gridfn(unknown_gfn=%d):\n" +" can't have want_ghost_zones && output_xyz_flag\n" +" && !radius_is_ghosted_flag!\n" +" unknown_radius_gfn=%d\n" +, + unknown_gfn, + unknown_radius_gfn); /*NOTREACHED*/ + +#ifdef CCTK_HDF5 + + using std::ostringstream; + using std::string; + using std::vector; + + herr_t herr; + + hid_t writer; + + // Get grid hierarchy extentsion from IOUtil + const ioGH * const iogh = (const ioGH *)CCTK_GHExtension (cctkGH, "IO"); + assert (iogh); + + struct stat fileinfo; + if (! iogh->recovered || stat(output_file_name, &fileinfo) != 0) { + writer = H5Fcreate (output_file_name, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + } else { + writer = H5Fopen (output_file_name, H5F_ACC_RDWR, H5P_DEFAULT); + } + assert (writer >= 0); + + hid_t datatype; + if (sizeof(fp) == sizeof(float)) { + datatype = H5T_NATIVE_FLOAT; + } else if (sizeof(fp) == sizeof(double)) { + datatype = H5T_NATIVE_DOUBLE; + } else if (sizeof(fp) == sizeof(long double)) { + datatype = H5T_NATIVE_LDOUBLE; + } else { + assert (0); + } + assert (datatype >= 0); + + for (int pn=0; pn<N_patches(); ++pn) { + const patch& p = ith_patch(pn); + + hsize_t shape[2]; + shape[0] = p.effective_N_isigma(want_ghost_zones); + shape[1] = p.effective_N_irho(want_ghost_zones); + + hid_t dataspace = H5Screate_simple (2, shape, NULL); + assert (dataspace >= 0); + + + + ostringstream datasetnamebuf; + datasetnamebuf << gfn_name + << " it=" << cctkGH->cctk_iteration + << " patch=" << pn; + string datasetnamestr = datasetnamebuf.str(); + assert (datasetnamestr.size() <= 256); // limit dataset name length + const char * const datasetname = datasetnamestr.c_str(); + assert (datasetname); + + ostringstream datasetnamebufx; + datasetnamebufx << "x" + << " it=" << cctkGH->cctk_iteration + << " patch=" << pn; + string datasetnamestrx = datasetnamebufx.str(); + assert (datasetnamestrx.size() <= 256); // limit dataset name length + const char * const datasetnamex = datasetnamestrx.c_str(); + assert (datasetnamex); + + ostringstream datasetnamebufy; + datasetnamebufy << "y" + << " it=" << cctkGH->cctk_iteration + << " patch=" << pn; + string datasetnamestry = datasetnamebufy.str(); + assert (datasetnamestry.size() <= 256); // limit dataset name length + const char * const datasetnamey = datasetnamestry.c_str(); + assert (datasetnamey); + + ostringstream datasetnamebufz; + datasetnamebufz << "z" + << " it=" << cctkGH->cctk_iteration + << " patch=" << pn; + string datasetnamestrz = datasetnamebufz.str(); + assert (datasetnamestrz.size() <= 256); // limit dataset name length + const char * const datasetnamez = datasetnamestrz.c_str(); + assert (datasetnamez); + + + + hid_t dataset = H5Dcreate (writer, datasetname, datatype, dataspace, H5P_DEFAULT); + assert (dataset >= 0); + + hid_t datasetx; + hid_t datasety; + hid_t datasetz; + if (output_xyz_flag) { + + datasetx = H5Dcreate (writer, datasetnamex, datatype, dataspace, H5P_DEFAULT); + assert (datasetx >= 0); + + datasety = H5Dcreate (writer, datasetnamey, datatype, dataspace, H5P_DEFAULT); + assert (datasety >= 0); + + datasetz = H5Dcreate (writer, datasetnamez, datatype, dataspace, H5P_DEFAULT); + assert (datasetz >= 0); + + } // if output_xyf_flag + + + + vector<fp> tmpdata; + vector<fp> tmpx, tmpy, tmpz; + tmpdata.resize (shape[0] * shape[1]); + if (output_xyz_flag) { + tmpx.resize (shape[0] * shape[1]); + tmpy.resize (shape[0] * shape[1]); + tmpz.resize (shape[0] * shape[1]); + } + + int elt = 0; + for (int irho = p.effective_min_irho(want_ghost_zones); + irho <= p.effective_max_irho(want_ghost_zones); + ++irho) + { + for (int isigma = p.effective_min_isigma(want_ghost_zones); + isigma <= p.effective_max_isigma(want_ghost_zones); + ++isigma) + { + + tmpdata.at(elt) = p.unknown_gridfn (ghosted_flag, unknown_gfn, irho,isigma); + if (output_xyz_flag) { + const fp r = p.unknown_gridfn(radius_is_ghosted_flag, + unknown_radius_gfn, + irho,isigma); + const fp rho = p.rho_of_irho(irho); + const fp sigma = p.sigma_of_isigma(isigma); + fp local_x, local_y, local_z; + p.xyz_of_r_rho_sigma(r,rho,sigma, + local_x,local_y,local_z); + const fp global_x = origin_x() + local_x; + const fp global_y = origin_y() + local_y; + const fp global_z = origin_z() + local_z; + tmpx.at(elt) = global_x; + tmpy.at(elt) = global_y; + tmpz.at(elt) = global_z; + } + ++elt; + + } + } + + + + herr = H5Dwrite (dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmpdata.front()); + assert (! herr); + + WriteAttribute (dataset, "name", gfn_name); + + herr = H5Dclose (dataset); + assert (! herr); + + + + if (output_xyz_flag) { + + herr = H5Dwrite (datasetx, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmpx.front()); + assert (! herr); + + WriteAttribute (datasetx, "name", "x"); + + herr = H5Dclose (datasetx); + assert (! herr); + + + + herr = H5Dwrite (datasety, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmpy.front()); + assert (! herr); + + WriteAttribute (datasety, "name", "y"); + + herr = H5Dclose (datasety); + assert (! herr); + + + + herr = H5Dwrite (datasetz, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tmpz.front()); + assert (! herr); + + WriteAttribute (datasetz, "name", "z"); + + herr = H5Dclose (datasetz); + assert (! herr); + + } // if output_xyz_flag + + + + herr = H5Sclose (dataspace); + assert (! herr); + + } // for pn + + + + herr = H5Fclose (writer); + assert (! herr); + +#else + error_exit(ERROR_EXIT, +"***** patch_system::output_unknown_gridfn(unknown_gfn=%d):\n" +" no HDF5 support compiled in. Cannot write output file \"%s\"\n!" +, + unknown_gfn, + output_file_name); /*NOTREACHED*/ + +#endif +} + +//****************************************************************************** //****************************************************************************** //****************************************************************************** |