aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/ChooseOutput.c12
-rw-r--r--src/Output.c2
-rw-r--r--src/Startup.c6
-rw-r--r--src/Write.c341
4 files changed, 274 insertions, 87 deletions
diff --git a/src/ChooseOutput.c b/src/ChooseOutput.c
index 7366fa7..04497c5 100644
--- a/src/ChooseOutput.c
+++ b/src/ChooseOutput.c
@@ -13,6 +13,7 @@
#include <string.h>
#include "cctk.h"
+#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "CactusBase/IOUtil/src/ioutil_Utils.h"
#include "ioJpegGH.h"
@@ -41,8 +42,6 @@ CCTK_FILEVERSION(CactusIO_IOJpeg_ChooseOutput_c)
/********************************************************************
******************** External Routines ************************
********************************************************************/
-void IOJpeg_ChooseOutput (const cGH *GH);
-
/*@@
@routine IOJpeg_ChooseOutput
@@ -60,13 +59,14 @@ void IOJpeg_ChooseOutput (const cGH *GH);
@vio in
@endvar
@@*/
-void IOJpeg_ChooseOutput (const cGH *GH)
+void IOJpeg_ChooseOutput (CCTK_ARGUMENTS)
{
int i, maxdim;
ioJpegGH *myGH;
int origin_index[3];
CCTK_REAL origin_phys[3];
- DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
GET_SLICE (out2D_xyplane_z, out_xyplane_z, origin_index[0], origin_phys[0]);
@@ -74,7 +74,7 @@ void IOJpeg_ChooseOutput (const cGH *GH)
GET_SLICE (out2D_yzplane_x, out_yzplane_x, origin_index[2], origin_phys[2]);
maxdim = CCTK_MaxDim ();
- myGH = (ioJpegGH *) CCTK_GHExtension (GH, "IOJpeg");
+ myGH = (ioJpegGH *) CCTK_GHExtension (cctkGH, "IOJpeg");
myGH->sp2xyz = (int **) malloc (3 * sizeof (int *));
for (i = 0; i < maxdim; i++)
@@ -83,7 +83,7 @@ void IOJpeg_ChooseOutput (const cGH *GH)
if (i > 0 && i < 3)
{
- IOUtil_2DPlanes (GH, i + 1, origin_index, origin_phys, myGH->sp2xyz[i]);
+ IOUtil_2DPlanes (cctkGH, i + 1, origin_index, origin_phys, myGH->sp2xyz[i]);
}
}
}
diff --git a/src/Output.c b/src/Output.c
index 312e4ca..ea72af8 100644
--- a/src/Output.c
+++ b/src/Output.c
@@ -377,7 +377,7 @@ static void SetOutputFlag (int vindex, const char *optstring, void *arg)
IOUtil_ParseOutputFrequency (CCTK_THORNSTRING, "IOJpeg::out_vars",
myGH->stop_on_parse_errors,
vindex, optstring,
- &myGH->out_every[vindex]);
+ &myGH->out_every[vindex], NULL);
}
}
}
diff --git a/src/Startup.c b/src/Startup.c
index 27566d8..c90472b 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -27,8 +27,6 @@ CCTK_FILEVERSION(CactusIO_IOJpeg_Startup_c)
/********************************************************************
******************** External Routines ************************
********************************************************************/
-void IOJpeg_Startup (void);
-
/********************************************************************
******************** Internal Routines ************************
@@ -49,13 +47,15 @@ static void *SetupGH (tFleshConfig *config, int conv_level, cGH *GH);
@calls CCTK_RegisterGHExtension
CCTK_RegisterGHExtensionSetupGH
@@*/
-void IOJpeg_Startup (void)
+int IOJpeg_Startup (void)
{
if (CCTK_MaxDim () >= 2)
{
CCTK_RegisterGHExtensionSetupGH (CCTK_RegisterGHExtension ("IOJpeg"),
SetupGH);
}
+
+ return 0;
}
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);
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-09 15:34:47 +0000