diff options
Diffstat (limited to 'src/Write.c')
| -rw-r--r-- | src/Write.c | 341 |
1 files changed, 264 insertions, 77 deletions
diff --git a/src/Write.c b/src/Write.c index 3fa5c7d..4ada69b 100644 --- a/src/Write.c +++ b/src/Write.c @@ -9,6 +9,9 @@ @@*/ +#include <assert.h> +#include <limits.h> +#include <math.h> #include <stdlib.h> #include <string.h> @@ -118,103 +121,287 @@ int IOJpeg_Write (const cGH *GH, CCTK_INT vindex, const char *alias) min = colormap_min; max = colormap_max; } - else + else if (CCTK_Equals (colormap, "auto")) { i = CCTK_LocalArrayReductionHandle ("minimum"); - i = CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_min); + CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_min); i = CCTK_LocalArrayReductionHandle ("maximum"); - i = CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_max); + CCTK_ReduceGridArrays (GH, 0, i, -1, 1, input_array , 1, input_array_type_codes, value_max); + } + else if (CCTK_Equals (colormap, "auto-old")) + { + i = CCTK_ReductionHandle ("minimum"); + CCTK_Reduce (GH, 0, i, 1, CCTK_VARIABLE_REAL, &min, 1, vindex); + i = CCTK_ReductionHandle ("maximum"); + CCTK_Reduce (GH, 0, i, 1, CCTK_VARIABLE_REAL, &max, 1, vindex); + } + else + { + CCTK_WARN (CCTK_WARN_ABORT, "Setting for parameter IOpjeg::colormap not recognised"); } /* get the number of slices to output */ /* in general: maxdir = gdata.dim * (gdata.dim - 1) / 2; */ maxdir = gdata.dim == 2 ? 1 : 3; - /* now do the actual I/O looping over all directions */ - for (dir = 0; dir < maxdir; dir++) + if (CCTK_EQUALS(gridpoints, "hyperslab")) { - /* get the directions to span the hyperslab */ - if (dir == 0) + + /* now do the actual I/O looping over all directions */ + for (dir = 0; dir < maxdir; dir++) { - dir_i = 0; dir_j = 1; /* xy */ - } - else if (dir == 1) - { - dir_i = 0; dir_j = 2; /* xz */ - } - else - { - dir_i = 1; dir_j = 2; /* yz */ - } + /* get the directions to span the hyperslab */ + if (dir == 0) + { + dir_i = 0; dir_j = 1; /* xy */ + } + else if (dir == 1) + { + dir_i = 0; dir_j = 2; /* xz */ + } + else + { + dir_i = 1; dir_j = 2; /* yz */ + } + + /* set the origin using the slice center from IOUtil */ + memset (origin, 0, sizeof (origin)); + if (gdata.grouptype == CCTK_GF) + { + origin[maxdir-dir-1] = myGH->sp2xyz[gdata.dim - 1][dir]; + } + + /* set the direction vector */ + memset (direction, 0, sizeof (direction)); + direction[dir_i] = direction[gdata.dim + dir_j] = 1; + + mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 2, + direction, origin, extent, + NULL, -1, NULL, hsize); + if (mapping < 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Failed to define hyperslab mapping for variable '%s'", + fullname); + continue; + } + { + // Sanity check + // (if this fails, we requested an insane number of grid points) + int max = INT_MAX; + int d; + for (d=0; d<2; ++d) { + assert (hsize[d] >= 0 && hsize[d] <= max); + if (hsize[d] > 0) max /= hsize[d]; + } + assert (CCTK_VarTypeSize (gdata.vartype) <= max); + } + total_hsize = hsize[0] * hsize[1] * CCTK_VarTypeSize (gdata.vartype); + assert (total_hsize >= 0); + if (total_hsize <= 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Selected hyperslab has zero size for variable '%s' " + "direction %d", fullname, dir); + Hyperslab_FreeMapping (mapping); + continue; + } + + /* allocate hyperslab buffer */ + hdata = myproc == 0 ? malloc (total_hsize) : NULL; + + /* get the hyperslab */ + i = Hyperslab_Get (GH, mapping, 0, vindex, 0, htype, hdata); + if (i) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Failed to extract hyperslab for variable '%s'", fullname); + } + + /* release the mapping structure */ + Hyperslab_FreeMapping (mapping); + + /* and dump the data to file */ + if (myproc == 0) + { + if (i == 0) + { + /* refine the 2D slice if requested */ + outdata = refinement_factor > 1 ? RefineData (hsize, hdata) : hdata; + if (! outdata) + { + CCTK_WARN (1, "Failed to refine 2D array"); + outdata = hdata; + } + WriteData (GH, vindex, alias, gdata.dim, dir, min, max, hsize, outdata); + if (outdata != hdata) + { + free (outdata); + } + } + + /* clean up */ + free (hdata); + } + } /* end of looping through xyz directions */ - /* set the origin using the slice center from IOUtil */ - memset (origin, 0, sizeof (origin)); - if (gdata.grouptype == CCTK_GF) + } + else if (CCTK_EQUALS(gridpoints, "interpolate")) + { + + int interpolator; + int options_table; + int coord_handle; + + CCTK_REAL * restrict coordsx; + CCTK_REAL * restrict coordsy; + CCTK_REAL * restrict coordsz; + CCTK_POINTER_TO_CONST coords[3]; + CCTK_INT * restrict inputs; + CCTK_INT * restrict output_types; + CCTK_POINTER * restrict outputs; + CCTK_INT * restrict operation_codes; + CCTK_INT * restrict time_deriv_order; + + const int nvars = 1; + int npoints; + + int n; + int i, j, k; + int d; + int ierr; + + CCTK_REAL *outdata; + + + + interpolator = CCTK_InterpHandle (interpolator_name); + assert (interpolator >= 0); + + options_table = Util_TableCreateFromString (interpolator_options); + assert (options_table >= 0); + + coord_handle = CCTK_CoordSystemHandle (interpolator_coordinates); + assert (coord_handle >= 0); + + + { - origin[maxdir-dir-1] = myGH->sp2xyz[gdata.dim - 1][dir]; + ierr = Util_TableSetInt (options_table, 1, "want_global_mode"); + assert (! ierr); } - - /* set the direction vector */ - memset (direction, 0, sizeof (direction)); - direction[dir_i] = direction[gdata.dim + dir_j] = 1; - - mapping = Hyperslab_GlobalMappingByIndex (GH, vindex, 2, - direction, origin, extent, - NULL, -1, NULL, hsize); - if (mapping < 0) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Failed to define hyperslab mapping for variable '%s'", - fullname); - continue; + + + + /* 2D Arrays */ + npoints = myproc==0 ? array2d_npoints_i * array2d_npoints_j : 0; + + coordsx = malloc (npoints * sizeof * coordsx); + assert (npoints==0 || coordsx); + coordsy = malloc (npoints * sizeof * coordsy); + assert (npoints==0 || coordsy); + coordsz = malloc (npoints * sizeof * coordsz); + assert (npoints==0 || coordsz); + coords[0] = coordsx; + coords[1] = coordsy; + coords[2] = coordsz; + + if (myproc==0) { + n = 0; + for (j=0; j<array2d_npoints_j; ++j) { + for (i=0; i<array2d_npoints_i; ++i) { + assert (n <= npoints); + coordsx[n] = array2d_x0 + i * array2d_dx_i + j * array2d_dx_j; + coordsy[n] = array2d_y0 + i * array2d_dy_i + j * array2d_dy_j; + coordsz[n] = array2d_z0 + i * array2d_dz_i + j * array2d_dz_j; + ++n; + } + } + assert (n == npoints); } - total_hsize = hsize[0] * hsize[1] * CCTK_VarTypeSize (gdata.vartype); - if (total_hsize <= 0) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Selected hyperslab has zero size for variable '%s' " - "direction %d", fullname, dir); - Hyperslab_FreeMapping (mapping); - continue; + + inputs = malloc (nvars * sizeof * inputs); + assert (inputs); + + for (n=0; n<nvars; ++n) { + inputs[n] = vindex; + if (inputs[n] < 0) { + inputs[n] = -1; + } } - - /* allocate hyperslab buffer */ - hdata = myproc == 0 ? malloc (total_hsize) : NULL; - - /* get the hyperslab */ - i = Hyperslab_Get (GH, mapping, 0, vindex, 0, htype, hdata); - if (i) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Failed to extract hyperslab for variable '%s'", fullname); + + output_types = malloc (nvars * sizeof * output_types); + assert (output_types); + + for (n=0; n<nvars; ++n) { + output_types[n] = CCTK_VARIABLE_REAL; } - - /* release the mapping structure */ - Hyperslab_FreeMapping (mapping); - - /* and dump the data to file */ - if (myproc == 0) - { - if (i == 0) - { - /* refine the 2D slice if requested */ - outdata = refinement_factor > 1 ? RefineData (hsize, hdata) : hdata; - if (! outdata) - { - CCTK_WARN (1, "Failed to refine 2D array"); - outdata = hdata; - } - WriteData (GH, vindex, alias, gdata.dim, dir, min, max, hsize, outdata); - if (outdata != hdata) - { - free (outdata); - } + + outputs = malloc (nvars * sizeof * outputs); + assert (outputs); + + outdata = malloc (npoints * sizeof * outdata); + assert (npoints==0 || outdata); + + for (n=0; n<nvars; ++n) { + outputs[n] = outdata; + } + + ierr = CCTK_InterpGridArrays + (GH, 3, interpolator, options_table, coord_handle, + npoints, CCTK_VARIABLE_REAL, coords, + nvars, inputs, + nvars, output_types, outputs); + assert (! ierr); + + + + /* Multiply values by r */ + if (myproc==0) { + for (n=0; n<npoints; ++n) { + const CCTK_REAL x = coordsx[n]; + const CCTK_REAL y = coordsy[n]; + const CCTK_REAL z = coordsz[n]; + const CCTK_REAL r = sqrt (pow (x, 2) + pow (y, 2) + pow (z, 2)); + outdata[n] *= r; } + } + + - /* clean up */ - free (hdata); + free (coordsx); + free (coordsy); + free (coordsz); + free (inputs); + free (output_types); + free (outputs); + + + + ierr = Util_TableDestroy (options_table); + assert (! ierr); + + + + if (myproc==0) { + const int dir = 0; + CCTK_INT hsize[2]; + hsize[0] = array2d_npoints_i; + hsize[1] = array2d_npoints_j; + + assert (CCTK_VarTypeI(vindex) == CCTK_VARIABLE_REAL); + + WriteData (GH, vindex, alias, gdata.dim, dir, + min, max, hsize, outdata); } - } /* end of looping through xyz directions */ + + free (outdata); + + } + else + { + CCTK_WARN (CCTK_WARN_ABORT, "unknown setting for parameter IOJpeg::gridpoints"); + } free (fullname); |