diff options
| -rw-r--r-- | par/jpeg_amr.par | 99 | ||||
| -rw-r--r-- | par/jpeg_amr.th | 32 | ||||
| -rw-r--r-- | param.ccl | 98 | ||||
| -rw-r--r-- | src/ChooseOutput.c | 12 | ||||
| -rw-r--r-- | src/Output.c | 2 | ||||
| -rw-r--r-- | src/Startup.c | 6 | ||||
| -rw-r--r-- | src/Write.c | 341 |
7 files changed, 498 insertions, 92 deletions
diff --git a/par/jpeg_amr.par b/par/jpeg_amr.par new file mode 100644 index 0000000..7d91c1b --- /dev/null +++ b/par/jpeg_amr.par @@ -0,0 +1,99 @@ +!DESC "Orbiting binary sources, example of jpeg IO with AMR" + +ActiveThorns = " + Boundary + CartGrid3D + CoordBase + IOUtil + InitBase + LocalInterp + SymBase + Time + + HTTPD + HTTPDExtra + Socket + + GSL + jpeg6b + + IOJpeg + + Carpet + CarpetIOBasic + CarpetInterp + CarpetLib + CarpetReduce + CarpetRegrid + CarpetSlab + LoopControl + + IDScalarWaveC + WaveBinarySource + WaveToyC + + AEILocalInterp +" + +Cactus::cctk_run_title = "WaveToy/httpd Example" +Cactus::terminate = "never" + +InitBase::initial_data_setup_method = "init_some_levels" +Carpet::init_fill_timelevels = yes + +Carpet::max_refinement_levels = 2 +CarpetRegrid::refinement_levels = 2 + +Carpet::domain_from_coordbase = yes +Grid::type = "coordbase" + +CoordBase::spacing = "numcells" +CoordBase::ncells_x = 70 +CoordBase::ncells_y = 70 +CoordBase::ncells_z = 70 +CoordBase::xmin = -1.0 +CoordBase::ymin = -1.0 +CoordBase::zmin = -1.0 +CoordBase::xmax = +1.0 +CoordBase::ymax = +1.0 +CoordBase::zmax = +1.0 + +IOJpeg::out_vars = "wavetoy::phi" +IOJpeg::mode = "remove" +IOJpeg::gridpoints = "interpolate" + +IOJpeg::out_every = 5 +IOJpeg::refinement_factor = 10 +#IOJpeg::colormap = "auto-old" # does not work with mesh refinement +IOJpeg::colormap = "custom" +IOJpeg::colormap_min = 0.0 # -1.0 +IOJpeg::colormap_max = 1.0 # +1.0 +IOJpeg::colormap_factor = 16 +IOJpeg::out_dir = $parfile + +IOJpeg::multiply_by_radius = yes +IOJpeg::array2d_x0 = -1.0 +IOJpeg::array2d_y0 = -1.0 +IOJpeg::array2d_z0 = 0.0 +IOJpeg::array2d_npoints_i = 101 +IOJpeg::array2d_dx_i = 0.02 +IOJpeg::array2d_dy_i = 0 +IOJpeg::array2d_dz_i = 0 +IOJpeg::array2d_npoints_j = 101 +IOJpeg::array2d_dx_j = 0 +IOJpeg::array2d_dy_j = 0.02 +IOJpeg::array2d_dz_j = 0 + +IOBasic::outInfo_every = 1 +IOBasic::outInfo_vars = "wavetoy::phi" + +Time::dtfac = 0.25 + +IDScalarWave::initial_data = "none" + +WaveBinarySource::binary_omega = 26 +WaveBinarySource::binary_charge = 0.0001 +WaveBinarySource::binary_radius = 0.25 +WaveBinarySource::binary_size = 0.1 + +WaveToy::bound = "radiation" diff --git a/par/jpeg_amr.th b/par/jpeg_amr.th new file mode 100644 index 0000000..336cfdc --- /dev/null +++ b/par/jpeg_amr.th @@ -0,0 +1,32 @@ +CactusBase/Boundary +CactusBase/CartGrid3D +CactusBase/CoordBase +CactusBase/IOUtil +CactusBase/InitBase +CactusBase/LocalInterp +CactusBase/SymBase +CactusBase/Time + +CactusConnect/HTTPD +CactusConnect/HTTPDExtra +CactusConnect/Socket + +CactusExternal/GSL +CactusExternal/jpeg6b + +CactusIO/IOJpeg + +Carpet/Carpet +Carpet/CarpetIOBasic +Carpet/CarpetInterp +Carpet/CarpetLib +Carpet/CarpetReduce +Carpet/CarpetRegrid +Carpet/CarpetSlab +Carpet/LoopControl + +CactusWave/IDScalarWaveC +CactusWave/WaveBinarySource +CactusWave/WaveToyC + +AEIThorns/AEILocalInterp @@ -37,6 +37,12 @@ KEYWORD mode "Output mode to use" STEERABLE = ALWAYS "standard" :: "Generate a file for each out_every timesteps" } "standard" +KEYWORD gridpoints "How to access grid points" STEERABLE = RECOVER +{ + "hyperslab" :: "use locations of grid points" + "interpolate" :: "interpolate to arbitrary points" +} "hyperslab" + ######################## # Specific to jpegs @@ -58,8 +64,9 @@ INT colormap_factor "How to scale float values to rgb color" STEERABLE = ALWAYS KEYWORD colormap "How to set the colormap" STEERABLE = ALWAYS { - "auto" :: "Set automatically using min/max of grid variables" - "custom" :: "Set min/max manually" + "auto" :: "Set automatically using min/max of grid variables" + "auto-old" :: "Set automatically using min/max of grid variables, using the old reduction interface which is still used by Carpet" + "custom" :: "Set min/max manually" } "custom" REAL colormap_min "minimum value to be mapped to colors" STEERABLE = ALWAYS @@ -77,9 +84,9 @@ INT refinement_factor "Refine each 2D slice by a certain factor (using interpola } 1 -################################ -# Choosing what planes to output -################################ +################################################### +# Choosing what planes to output when hyperslabbing +################################################### REAL out2D_yzplane_x "x-coord for 2D planes in yz" STEERABLE = RECOVER { *:* :: "A value between [xmin, xmax]" @@ -113,6 +120,87 @@ INT out2D_xyplane_zi "z-index (from 0) for 2D planes in xy" STEERABLE = RECOVER } -1 +################################################### +# Choosing what region to output when interpolating +################################################### + +STRING interpolator_name "Name of the interpolator" STEERABLE=always +{ + ".*" :: "must be a registered interpolator" +} "Lagrange polynomial interpolation" + +STRING interpolator_options "Options for the interpolator" STEERABLE=always +{ + ".*" :: "must be a valid option specification" +} "order=2" + +STRING interpolator_coordinates "Coordinate system" STEERABLE=always +{ + ".*" :: "must be a registered coordinate system" +} "cart3d" + +BOOLEAN multiply_by_radius "Multiply valus by r" STEERABLE=always +{ +} "no" + + + +REAL array2d_x0 "Origin" STEERABLE=always +{ + *:* :: "" +} 0.0 + +REAL array2d_y0 "Origin" STEERABLE=always +{ + *:* :: "" +} 0.0 + +REAL array2d_z0 "Origin" STEERABLE=always +{ + *:* :: "" +} 0.0 + +INT array2d_npoints_i "Number of grid points for the 2D grid arrays in the i direction" +{ + 0:* :: "" +} 10 + +REAL array2d_dx_i "Spacing" STEERABLE=always +{ + 0.0:* :: "" +} 0.0 + +REAL array2d_dy_i "Spacing" STEERABLE=always +{ + 0.0:* :: "" +} 0.0 + +REAL array2d_dz_i "Spacing" STEERABLE=always +{ + 0.0:* :: "" +} 0.0 + +INT array2d_npoints_j "Number of grid points for the 2D grid arrays in the j direction" +{ + 0:* :: "" +} 10 + +REAL array2d_dx_j "Spacing" STEERABLE=always +{ + 0.0:* :: "" +} 0.0 + +REAL array2d_dy_j "Spacing" STEERABLE=always +{ + 0.0:* :: "" +} 0.0 + +REAL array2d_dz_j "Spacing" STEERABLE=always +{ + 0.0:* :: "" +} 0.0 + + ############################################################################# ### import IOUtil parameters ############################################################################# 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); |