diff --git a/src/rz_array.c b/src/rz_array.c new file mode 100644 index 0000000000000000000000000000000000000000..828bfc3768d347e692a30724480823f9399cb1ee --- /dev/null +++ b/src/rz_array.c @@ -0,0 +1,444 @@ +/** @file rz_array.c + * @brief General functions to work with FORTRAN-compatible 2d/3d arrays. + */ +#include +#include +#include + +#include "parameters.h" +#include "proto.h" +#include "rz_array.h" +#include "species.h" + +long int used_disk_space = 0; +/*!< This is the total disk space used so far for the array dumps. + It is used to limit the disk space used by a simulation and thus + avoid that one simulation gets crazy and stops other simulations + that are sharing the same disk space. */ + +static void check_disk_space(void); + +/** @brief Creates a new 3D array structure. */ +rz_array_t * +rz_new_3d_a (int r0, int z0, int r1, int z1, int ntheta) +{ + int rmax, zmax; + rz_array_t *array; + + debug (3, "rz_new_3d_a (%d, %d, %d, %d, %d)\n", r0, z0, r1, z1, ntheta); + + /* We allow a margin of two cells to set the boundary conditions (see below)*/ + r0 -= 2; z0 -= 2; r1 += 2; z1 += 2; + + rmax = r1 - r0; + zmax = z1 - z0; + + assert (rmax > 0 && zmax > 0 && ntheta > 0); + + debug (3, "r0 = %d, z0 = %d, r1 = %d, z1 = %d\n", r0, z0, r1, z1); + + array = (rz_array_t *) xmalloc (sizeof(rz_array_t)); + + if (ntheta != 1) { + array->len = rmax * zmax * (ntheta + 4); + array->theta0 = -2; + } else { + array->len = rmax * zmax; + array->theta0 = 0; + } + + /* Initially all arrays are zeroed. */ + array->data = (REAL *) xcalloc (array->len, sizeof(REAL)); + + + /* These strides are to make the array compatible with array(ir, iz) + in FORTRAN. */ + array->strides[R_INDX] = 1; + array->strides[Z_INDX] = rmax; + array->strides[THETA_INDX] = rmax * zmax; + + array->dim = ntheta == 1? 2: 3; + + debug (3, "rmax = %d\n", rmax); + debug (3, "strides = {%d, %d}\n", array->strides[0], array->strides[1]); + + array->r0 = r0; + array->z0 = z0; + /* Note that these are NOT the r0, z0 we received. The idea here is that + RZ(array, r0, z0) will produce the correct result, where the allocated + array looks like (in the picture, r0, z0 ARE the received ones). + Note: in the picture there is only one buffer cell, while we actually + allocate two. + + +--------------+--------------+...+--------------+--------------+ + |r0 - 1, z0 - 1| | | | | + +--------------+--------------+...+--------------+--------------+ + | | r0, z0 | | | | + +--------------+--------------+...+--------------+--------------+ + ... |...| ... + +--------------+--------------+...+--------------+--------------+ + | | | |r1 - 1, z1 - 1| | + +--------------+--------------+...+--------------+--------------+ + | | | | | r1, z1 | + +--------------+--------------+...+--------------+--------------+ + + but the rows z0 +/- 1 and the columns r0 +/- 1 do not belong to the + physical space and are used only to specify boundary conditions. + */ + + array->nr = r1 - r0 - 4; + array->nz = z1 - z0 - 4; + array->ntheta = ntheta; + + array->host = NULL; + + return array; +} + +/** @brief Creates guest array: This is an array that contains + * part of another one. + * + * We use the memory of the original array and do not allocate + * new memory. + */ +rz_array_t * +rz_guest (rz_array_t *host, int r0, int z0, int r1, int z1) +{ + int ntheta, rmax, zmax; + rz_array_t *array; + + ntheta = host->ntheta; + array = (rz_array_t *) xmalloc (sizeof(rz_array_t)); + + rmax = r1 - r0; + zmax = z1 - z0; + + if (ntheta != 1) { + array->len = rmax * zmax * (ntheta + 4); + } else { + array->len = rmax * zmax; + } + + array->r0 = r0; + array->z0 = z0; + + array->theta0 = host->theta0; + + array->strides[R_INDX] = host->strides[R_INDX]; + array->strides[Z_INDX] = host->strides[Z_INDX]; + array->strides[THETA_INDX] = host->strides[THETA_INDX]; + + array->dim = host->dim; + + array->nr = r1 - r0 - 4; + array->nz = z1 - z0 - 4; + array->ntheta = host->ntheta; + + array->data = RZTP (host, r0, z0, array->theta0); + array->host = host; + + return array; +} + +/** @brief Sets the data of an array to zero */ +void +rz_set_zero (rz_array_t *array) +{ + debug (3, "rz_set_zero\n"); + memset (array->data, 0, array->len * sizeof(REAL)); +} + +/** @brief Sets periodic boundary conditions for theta in an array.*/ +void +rz_set_periodic (rz_array_t *array) +{ + /* If the array is 2D, we do nothing. */ + if (array->ntheta == 1) return; + + rz_copy_modes (array, array->r0, array->z0, + array, array->r0, array->z0, + array->nr + 4, array->nz + 4, array->ntheta - 1, -1); + + rz_copy_modes (array, array->r0, array->z0, + array, array->r0, array->z0, + array->nr + 4, array->nz + 4, array->ntheta - 2, -2); + + rz_copy_modes (array, array->r0, array->z0, + array, array->r0, array->z0, + array->nr + 4, array->nz + 4, 0, array->ntheta); + + rz_copy_modes (array, array->r0, array->z0, + array, array->r0, array->z0, + array->nr + 4, array->nz + 4, 1, array->ntheta + 1); + +} + + +/** @brief Sets boundaries for one array reading data from other array + * (but both can be the same: see rz_set_bnd). + * + * Starting at *start_xxx, it sweeps a 2-dimensional subspace of the grid + * in dimensions dim1 and dim2 (and perpendicular to dim0). There + * it sets the boudary conditions. + * + * The value of from is multiplied by sign. + * Thus, if the boundary is itself, is -1 for homogeneous dirichlet + * and 1 for homogeneous Neumann (though you should better use + * BND_CND_HNEUMANN and BND_CND_HDIRICHLET) + * + * @a inout is -1 if we look/set the values on smaller values + * of dim0 1 if we look at larger values. + * (though it is better to use BND_INWARD=-1 and BND_OUTWARD=1). + * + * Note: it is better to call this function _before_ rz_set_periodic. + */ +void +rz_copy_bnd (rz_array_t *from, rz_array_t *to, + int sign, REAL *start_from, REAL *start_to, + int dim0, int inout_from, int inout_to, + int dim1, int dim1_0, int dim1_1, + int dim2, int dim2_0, int dim2_1) +{ + int i, j; + REAL *pfrom, *pto; + + for (i = dim1_0; i < dim1_1; i++) { + for (j = dim2_0; j < dim2_1; j++) { + pfrom = start_from + i * from->strides[dim1] + j * from->strides[dim2]; + pto = start_to + i * to->strides[dim1] + j * to->strides[dim2]; + + *(pto + inout_to * to->strides[dim0]) = sign * (*pfrom); + + *(pto + 2 * inout_to * to->strides[dim0]) = + sign * (*(pfrom + inout_from * from->strides[dim0])); + } + } +} + + +/** @brief Sets boundary conditions (Neumann/Dirichlet) reading from + * the array itself. + * + * For example, to set Neumann conditions at r = 0 (dim0 = R_INDX), + * we call it as + * + * rz_set_bnd (array, 1, RZTP (array, 0, 0, 0), R_INDX, 1, + * Z_INDX, z0, z1, THETA_INDX, 0, ntheta); + */ +void +rz_set_bnd (rz_array_t *array, int sign, REAL *astart, int dim0, int inout, + int dim1, int dim1_0, int dim1_1, + int dim2, int dim2_0, int dim2_1) +{ + rz_copy_bnd (array, array, sign, astart, astart, dim0, -inout, inout, + dim1, dim1_0, dim1_1, dim2, dim2_0, dim2_1); +} + + +/** @brief Frees an array */ +void +rz_free (rz_array_t *array) +{ + debug (3, "rz_free\n"); + + if (NULL == array->host) free (array->data); + free (array); +} + +/** @brief Copies a sub-array of @a rn x @a zn from @a fro + * (starting at (@a rfro, z@a fro)) to @a to (starting at + * (@a rto, @a zto)). */ +void +rz_copy (rz_array_t *fro, int rfro, int zfro, + rz_array_t *to, int rto, int zto, + int rn, int zn) +{ + int i; + REAL *pfro, *pto; + + debug (3, "rz_copy\n"); + + pfro = RZP (fro, rfro, zfro); + pto = RZP (to, rto, zto); + + for (i = 0; i < zn; i++) { + memcpy (pto, pfro, sizeof(REAL) * rn); + + pfro += fro->strides[Z_INDX]; + pto += to->strides[Z_INDX]; + } +} + +/** @brief Copies a sub-array of the mode nmode of @a rn x @a zn from + * @a fro (starting at (@a rfro, @a zfro)) to @a to (starting at + * (@a rto, @a zto)). +*/ +void +rz_copy_modes (rz_array_t *fro, int rfro, int zfro, + rz_array_t *to, int rto, int zto, + int rn, int zn, int nmode_fro, int nmode_to) +{ + int i, j; + REAL *pfro, *pto; + + debug (3, "rz_copy\n"); + + pfro = RZTP (fro, rfro, zfro, nmode_fro); + pto = RZTP (to, rto, zto, nmode_to); + + for (i = 0; i < zn; i++) { + /* memcpy (pto, pfro, sizeof(REAL) * rn); Maybe not thread safe? */ + for (j = 0; j < rn; j++ ) *(pto + j) = *(pfro + j); + + pfro += fro->strides[Z_INDX]; + pto += to->strides[Z_INDX]; + } +} + + +/** @brief Reads or writes a @a rz_array, depending on the string mode, + * "r" for reading, "w" for writing. + * + * Note however that we must already have the dimensions of the array in it. + */ +void +rz_dump (rz_array_t *rz_array, const char *fname, const char *mode, + int r0, int z0, int r1, int z1) +{ + FILE *fp; + int ir, iz; + int writing; + double f; + + debug (3, "rz_dump(\"%s\", \"%s\", %d, %d, %d, %d)\n", fname, mode, + r0, z0, r1, z1); + + if (0 == strcmp (mode, "w")) { + writing = TRUE; + } else if (0 == strcmp (mode, "r")) { + writing = FALSE; + } else { + fatal ("Unknown mode for rz_dump\n"); + } + + fp = fopen (fname, mode); + + if (fp == NULL) { + warning("Unable to open %s for %s\n", fname, + writing? "writing": "reading"); + return; + } + + for (ir = r0; ir < r1; ir++) + for (iz = z0; iz < z1; iz++) { + if (writing) + fprintf (fp, "%15.5e\n", RZ(rz_array, ir, iz)); + else { + if (fscanf (fp, "%lf", &f) != 1) { + warning ("Error reading file %s, at line %d, (ir, iz) = (%d, %d)\n", + fname, iz + ir * (z1 - z0), ir, iz); + f = 0.0; + } + *RZP (rz_array, ir, iz) = f; + } + } + + used_disk_space += ftell (fp); + fclose(fp); + + check_disk_space (); +} + +/** @brief rz_dump_3d QQQQ */ +void +rz_dump_3d (rz_array_t *rz_array, const char *fname, const char *mode, + int r0, int z0, int r1, int z1, int ntheta) +{ + FILE *fp; + int ir, iz, itheta; + int writing; + double f; + + debug (3, "rz_dump_3d(\"%s\", \"%s\", %d, %d, %d, %d, %d)\n", + fname, mode, r0, z0, r1, z1, ntheta); + + if (0 == strcmp (mode, "w")) { + writing = TRUE; + } else if (0 == strcmp (mode, "r")) { + writing = FALSE; + } else { + fatal ("Unknown mode for rz_dump_3d\n"); + } + + fp = fopen (fname, mode); + + if (fp == NULL) { + warning ("Unable to open %s for %s\n", fname, + writing? "writing": "reading"); + return; + } + + for (itheta = 0; itheta < ntheta; itheta++) { + for (ir = r0; ir < r1; ir++) { + for (iz = z0; iz < z1; iz++) { + if (writing) + fprintf (fp, "%15.5e\n", RZT (rz_array, ir, iz, itheta)); + else { + if (fscanf (fp, "%lf", &f) != 1) { + warning ("Error reading file %s, at line %d, " + "(ir, iz, itheta) = (%d, %d, %d)\n", + fname, iz + (ir + itheta * (r1 - r0)) * (z1 - z0), + ir, iz, itheta); + f = 0.0; + } + *RZTP (rz_array, ir, iz, itheta) = f; + } + } + } + } + + /* Add the disk space used to the total disk space that we have used so far + */ + + used_disk_space += ftell (fp); + fclose(fp); + + check_disk_space (); +} + +/** @brief rz_axis_dump QQQQ */ +void + rz_axis_dump (const char *fname, int x0, int x1, double delta) +{ + FILE *fp; + int i; + + debug (3, "rz_axis_dump(\"%s\", %d, %d, %f)\n", fname, x0, x1, delta); + + fp = fopen (fname, "w"); + + if (fp == NULL) { + warning ("Unable to open %s\n", fname); + return; + } + + for (i = x0; i < x1; i++) + fprintf (fp, "%15.5e\n", ((double) i + 0.5) * delta); + + + used_disk_space += ftell (fp); + fclose(fp); + + check_disk_space (); +} + +/** @brief check_disk_space QQQQ */ +static void +check_disk_space (void) +{ + if (used_disk_space > (((long int) max_disk_space_mb) << 20)) { + fatal ("The disk space limit has been surpassed. " + "Increase max_disk_space_mb if you really need more space\n"); + } +} +