CHarm#
This chapter demonstrates how to perform some basic tasks with
CHarm. For a detailed description of the functions to be used, see the
CHarm (C API) chapter. The source codes of the examples can be found in the
cookbook/c directory.
A single header file of the entire library is available as:
charm/charmf.hfor single precision,charm/charm.hfor double precision, andcharm/charmq.hfor quadruple precision.
You may combine various precisions in a single program, provided that you follow the rules from the CHarm in single and quadruple precision chapter below. It is not recommended though to combine various releases (e.g., single precision from v0.0.0 and double precision form v0.1.0).
The chapters that follow show how to work with CHarm in double, single and quadruple precision, respectively.
CHarm in double precision#
In this section, we assume you have successfully installed CHarm (see Installing) in double precision with OpenMP enabled into the default installation paths, and that you know how to compile a C code and link external libraries. A brief example of the compilation under Linux with GCC is provided in the Compilation on Linux section below.
Compilation on Linux#
Assuming your working directory is cookbook/c, an example compilation
with the static library might look like this:
gcc -fopenmp shcs.c -l:libcharm.a -lfftw3 -lfftw3_omp -lm
Compilation using the shared library:
gcc -fopenmp -Wl,-rpath -Wl,/usr/local/lib shcs.c \
-lcharm -lfftw3 -lfftw3_omp -lm
After a successful compilation, you are ready to execute the compiled binary as
./a.out
You should see the following output:
C( 9, 0) = 2.7671430085300001e-08
S( 9, 0) = 0.0000000000000000e+00
C( 9, 4) = -9.0017922533600003e-09
S( 9, 4) = 1.9466677947499999e-08
C( 9, 9) = -4.7747538613200003e-08
S( 9, 9) = 9.6641284771399995e-08
Degree variance for harmonic degree 0 = 1.0000000000000000e+00
Degree variance for harmonic degree 4 = 2.5183167531865940e-12
Degree variance for harmonic degree 10 = 1.2631494966213168e-13
Difference degree variance for harmonic degree 0 = 0.0000000000000000e+00
Difference degree variance for harmonic degree 4 = 0.0000000000000000e+00
Difference degree variance for harmonic degree 10 = 0.0000000000000000e+00
Great, all done!
Note
Various approaches exist to compile the code using static/shared libraries. Use whatever works best for you.
Working with spherical harmonic coefficients#
This example shows how to read/write spherical harmonic coefficients from/to files and how to compute (difference) degree variances.
#include <stdio.h>
#include <stdlib.h>
#include <charm/charm.h>
int main(void)
{
/* INPUTS */
/* ===================================================================== */
/* Define the path to an input "gfc" text file with spherical harmonic
* coefficients. For details on the structure of the "gfc" file, see the
* description of the "charm_shc_read_gfc" function in the "charm_shc"
* module. */
char shcs_in_file[] = "../../data/input/EGM96-degree10.gfc";
/* Maximum harmonic degree to initialize, read and write spherical harmonic
* coefficients */
unsigned long nmax = 10;
/* Define the path to an output binary file with spherical harmonic
* coefficients. For details on the structure of the output binary file,
* see the description of the "charm_shc_write_bin" function in the
* "charm_shc" module. */
char shcs_out_file[] = "../../data/output/EGM96-degree10.shcs";
/* ===================================================================== */
/* ===================================================================== */
/* Read the spherical harmonic coefficients */
/* --------------------------------------------------------------------- */
/* Initialize a "charm_shc" structure up to degree "nmax". All
* coefficients will be initialized to zero and both the scaling constant
* "shcs->mu" and the radius of the reference sphere "shcs->r" will be set
* to "1.0". These values will later be overwritten when reading the input
* file. In case of failure, returned is "NULL", so always check the
* returned pointer for the "NULL" value! */
charm_shc *shcs = charm_shc_calloc(nmax, 1.0, 1.0);
if (shcs == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure");
exit(CHARM_FAILURE);
}
/* Many CHarm functions take a "charm_err" structure as their last
* parameter. This is to allow the called function to report an error (if
* encountered) back to the caller, including some useful information, such
* as the error message, error code, function name and so on. The error
* structure is defined in the "charm_err" module and can be initialized
* like this: */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* Now let's read spherical harmonic coefficients, the scaling constant and
* the radius of the reference sphere from the input "gfc" file. We are
* reading a static "gfc" model, so the "epoch" parameter is set to "NULL".
* */
charm_shc_read_gfc(shcs_in_file, nmax, NULL, shcs, err);
/* At this point, "shcs" should store the loaded spherical harmonic
* coefficients. However, before we continue in our program, we should
* check whether the previous function call was successful or not. To this
* end, we shall look into the "err" structure by calling an error handler.
* The first parameter of the error handler is the error structure itself,
* "err". The other parameter says that if there is indeed an error
* message in "err", the program prints details on the error and
* subsequently terminates. If you want to print the error but do not want
* to terminate your program, replace "1" by "0" (further details in the
* "charm_err" module). After you call the error handler, you can reuse
* the same "err" structure in your program, since the function resets the
* error structure to the default empty values.
*
* It is not absolutely necessary to call the error handler, but it is
* really really really recommended to always do so. */
charm_err_handler(err, 1);
/* --------------------------------------------------------------------- */
/* Now print some more or less randomly chosen spherical harmonic
* coefficients */
/* --------------------------------------------------------------------- */
/* Let's start with zonal coefficients */
unsigned long n = 9; /* Harmonic degree */
unsigned long m = 0; /* Harmonic order */
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs->s[m][n - m]);
/* Now some tesseral coefficients */
n = 9;
m = 4;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs->s[m][n - m]);
/* And finally some sectorial coefficients */
n = 9;
m = 9;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs->s[m][n - m]);
/* --------------------------------------------------------------------- */
/* Now let's save the coefficients to a binary file and then read them back
* to another structure for spherical harmonic cofficients. Just for the
* fun... */
/* --------------------------------------------------------------------- */
/* Write the coefficients. */
charm_shc_write_bin(shcs, nmax, shcs_out_file, err);
/* Again, we should call the error handler to see whether the previous call
* was successful or not. */
charm_err_handler(err, 1);
/* And now read back the coefficients to a new "charm_shc" structure called
* "shcs2" */
charm_shc *shcs2 = charm_shc_calloc(nmax, 1.0, 1.0);
if (shcs2 == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.");
exit(CHARM_FAILURE);
}
charm_shc_read_bin(shcs_out_file, nmax, shcs2, err);
charm_err_handler(err, 1);
/* --------------------------------------------------------------------- */
/* Compute degree variances from "shcs" */
/* --------------------------------------------------------------------- */
/* Compute degree variances of the input signal */
double *dv = (double *)malloc((nmax + 1) * sizeof(double));
if (dv == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
charm_shc_dv(shcs, nmax, dv, err);
charm_err_handler(err, 1);
/* Print some degree variances */
n = 0;
printf("\n\n");
printf("Degree variance for harmonic degree %lu = %0.16e\n", n, dv[n]);
n = 4;
printf("Degree variance for harmonic degree %lu = %0.16e\n", n, dv[n]);
n = 10;
printf("Degree variance for harmonic degree %lu = %0.16e\n", n, dv[n]);
/* --------------------------------------------------------------------- */
/* Now check whether "shcs" and "shcs2" contain the same coefficients by
* computing difference degree variances */
/* --------------------------------------------------------------------- */
double *ddv = (double *)malloc((nmax + 1) * sizeof(double));
if (ddv == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
charm_shc_ddv(shcs, shcs2, nmax, ddv, err);
charm_err_handler(err, 1);
/* Print some difference degree variances */
n = 0;
printf("\n\n");
printf("Difference degree variance for harmonic degree %lu = %0.16e\n", n,
ddv[n]);
n = 4;
printf("Difference degree variance for harmonic degree %lu = %0.16e\n", n,
ddv[n]);
n = 10;
printf("Difference degree variance for harmonic degree %lu = %0.16e\n", n,
ddv[n]);
/* --------------------------------------------------------------------- */
/* --------------------------------------------------------------------- */
/* Remember: always free the memory associated with the CHarm structures
* via the special CHarm "charm_*_free*" functions. Calling the usual
* "free" function will not deallocate the memory properly and will lead to
* memory leaks. */
charm_err_free(err);
charm_shc_free(shcs);
charm_shc_free(shcs2);
free(dv), free(ddv);
/* --------------------------------------------------------------------- */
printf("\nGreat, all done!\n");
exit(CHARM_SUCCESS);
/* ===================================================================== */
}
Spherical harmonic synthesis and analysis#
This section shows several examples.
A closed-loop test of analysis and synthesis with point data values. First, spherical harmonic coefficients are loaded from a text file. Then, they are used to synthesize a signal, from which a new set of coefficients is finally computed by surface spherical harmonic analysis. The two coefficients sets are then finally compared (should be equal, that is, the difference degree amplitudes should be negligibly small).
A solid spherical harmonic synthesis of point values at scattered points.
A solid spherical harmonic synthesis of point values at a custom grid of points.
A solid spherical harmonic synthesis of block-mean values at scattered cells.
A solid spherical harmonic synthesis of block mean values at a grid of cells.
Surface spherical harmonic analysis with block-mean values in cells.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <charm/charm.h>
int main(void)
{
/* INPUTS */
/* ===================================================================== */
/* Define the path to an input text file with spherical harmonic
* coefficients. For details on the structure of the text file, see the
* description of the "charm_shc_read_mtx" function in the "charm_shc"
* module. */
char shcs_in_file[] = "../../data/input/EGM96-degree10-mtx.txt";
/* Maximum harmonic degree of coefficients to read. The same degree is
* used later in the harmonic synthesis and harmonic analysis. */
unsigned long nmax = 10;
/* ===================================================================== */
/* ===================================================================== */
printf("Closed-loop experiment -- Spherical harmonic synthesis and "
"analysis\n");
printf("===========================\n");
/* Initialize a "charm_shc" structure */
charm_shc *shcs = charm_shc_calloc(nmax, 1.0, 1.0);
if (shcs == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.\n");
exit(CHARM_FAILURE);
}
/* Initialize a "charm_err" structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* Read spherical harmonic coefficients from the input text file. */
charm_shc_read_mtx(shcs_in_file, nmax, shcs, err);
charm_err_handler(err, 1);
/* Now let's say we do not want to use the zero-degree term. This can be
* achieved simply by setting the respective "C00" coefficient to zero. */
shcs->c[0][0] = 0.0;
/* Compute the Gauss--Legendre grid for a given "nmax" on a sphere that is
* "1000.0" metres above the reference sphere of spherical harmonic
* coefficients. */
charm_point *grd_pnt = charm_crd_point_gl(nmax, shcs->r + 1000.0);
if (grd_pnt == NULL)
{
fprintf(stderr, "Failed to compute the Gauss--Legendre grid.\n");
exit(CHARM_FAILURE);
}
/* Allocate the memory for the synthesized signal */
double *f = (double *)malloc(grd_pnt->npoint * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Perform the synthesis. Since the Gauss--Legendre grid resides "1000.0"
* metres above the reference sphere of the coefficients, performed is
* solid spherical harmonic synthesis. */
charm_shs_point(grd_pnt, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Print some synthesized values */
printf("Print some synthesized values of the signal...\n");
size_t i = 0;
size_t j = 0;
printf("f(%zu, %zu) = %0.16e\n", i, j, f[i * grd_pnt->nlon + j]);
i = grd_pnt->nlat / 2;
j = grd_pnt->nlon / 2;
printf("f(%zu, %zu) = %0.16e\n", i, j, f[i * grd_pnt->nlon + j]);
/* Initialize a new structure of spherical harmonic coefficients for
* coefficients to be computed by the harmonic analysis. Used are the same
* scalling constants as in "shcs". */
charm_shc *shcs2 = charm_shc_calloc(nmax, shcs->mu, shcs->r);
if (shcs2 == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.\n");
exit(CHARM_FAILURE);
}
/* Now use the synthesized signal and compute back its spherical harmonic
* coefficients by harmonic analysis. The output coefficients in "shcs2"
* should be the same as the input ones "shcs". Note that the signal "f"
* resides on a sphere that is "1000.0" metres above the reference sphere.
* Therefore, the coefficients are automatically properly rescaled to the
* desired sphere defined by "shcs2->r". */
charm_sha_point(grd_pnt, f, nmax, shcs2, err);
charm_err_handler(err, 1);
/* Now check whether "shcs" and "shcs2" are the same by computing their
* difference degree amplitudes */
double *dda = (double *)malloc((nmax + 1) * sizeof(double));
if (dda == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
charm_shc_dda(shcs, shcs2, nmax, dda, err);
charm_err_handler(err, 1);
/* Print some difference degree amplitudes */
unsigned long n = 0;
printf("\n\n");
printf("Now print the difference degree amplitudes. ");
printf("These should be very small, say, at the order of 1e-18 "
"or less...\n");
printf("Difference degree amplitude for harmonic degree %lu = %0.16e\n",
n, dda[n]);
n = 4;
printf("Difference degree amplitude for harmonic degree %lu = %0.16e\n",
n, dda[n]);
n = 10;
printf("Difference degree amplitude for harmonic degree %lu = %0.16e\n",
n, dda[n]);
/* We will not need some of the structures and arrays anymore, so let's
* free them. */
charm_crd_point_free(grd_pnt);
charm_shc_free(shcs2);
free(f);
free(dda);
printf("===========================\n\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Solid spherical harmonic synthesis at scattered points\n");
printf("===========================\n");
/* Let's create a "charm_point" structure to store 3 scattered points. */
charm_point *sctr_pnt = charm_crd_point_calloc(CHARM_CRD_POINT_SCATTERED,
3, 3);
if (sctr_pnt == NULL)
{
fprintf(stderr,
"Failed to initialize the \"charm_point\" structure.\n");
exit(CHARM_FAILURE);
}
/* Now let's feed "sctr_pnt" with some computation points (angular values
* must be provided in radians). */
sctr_pnt->lat[0] = 0.1;
sctr_pnt->lat[1] = 0.436231;
sctr_pnt->lat[2] = -0.9651;
sctr_pnt->lon[0] = 0.0;
sctr_pnt->lon[1] = 1.53434;
sctr_pnt->lon[2] = 4.2316;
sctr_pnt->r[0] = shcs->r + 1000.0;
sctr_pnt->r[1] = shcs->r + 2000.0;
sctr_pnt->r[2] = shcs->r + 3000.0;
/* Allocate memory to store the synthesized signal */
f = (double *)malloc(sctr_pnt->npoint * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_point(sctr_pnt, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* It's really this easy! */
/* Print the synthesized values */
for (size_t i = 0; i < sctr_pnt->npoint; i++)
printf("%0.16e\n", f[i]);
charm_crd_point_free(sctr_pnt);
free(f);
printf("===========================\n\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Solid spherical harmonic synthesis at a custom grid of points\n");
printf("===========================\n");
/* Initialize a charm_point structure to hold a "5 x 10" point grid. */
grd_pnt = charm_crd_point_calloc(CHARM_CRD_POINT_GRID, 5, 10);
if (grd_pnt == NULL)
{
fprintf(stderr,
"Failed to initialize the \"charm_point\" structure.\n");
exit(CHARM_FAILURE);
}
/* Define some grid points */
for (size_t i = 0; i < grd_pnt->nlat; i++)
{
grd_pnt->lat[i] = M_PI_2 - ((double)i / (double)grd_pnt->nlat) * M_PI;
grd_pnt->r[i] = shcs->r + (double)i;
}
for (size_t j = 0; j < grd_pnt->nlon; j++)
{
grd_pnt->lon[j] = ((double)j / (double)grd_pnt->nlon) * (2.0 * M_PI);
}
/* Initialize an array to store the synthesized signal */
f = (double *)malloc(grd_pnt->npoint * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_point(grd_pnt, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Print the synthesized values */
for (size_t i = 0; i < grd_pnt->nlat; i++)
for (size_t j = 0; j < grd_pnt->nlon; j++)
printf("%0.16e\n", f[i * grd_pnt->nlon + j]);
charm_crd_point_free(grd_pnt);
free(f);
printf("===========================\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Solid spherical harmonic synthesis at scattered cells\n");
printf("===========================\n");
/* Initialize a "charm_cell" structure to hold 3 scattered cells. */
charm_cell *sctr_cell = charm_crd_cell_calloc(CHARM_CRD_CELL_SCATTERED,
3, 3);
if (sctr_cell == NULL)
{
fprintf(stderr, "Failed to initialize scattered cells.\n");
exit(CHARM_FAILURE);
}
/* Define more or less randomly some scattered cells */
/* --------------------------------------------------------------------- */
/* The maximum latitude of the first cell */
sctr_cell->latmax[0] = 0.323413435;
/* The minimum latitude of the first cell */
sctr_cell->latmin[0] = sctr_cell->latmax[0] - 0.234323;
/* The maximum latitude of the second cell */
sctr_cell->latmax[1] = -0.90234320952;
/* The minimum latitude of the second cell */
sctr_cell->latmin[1] = sctr_cell->latmax[1] - 0.4456;
/* And so on... */
sctr_cell->latmax[2] = 0.0;
sctr_cell->latmin[2] = sctr_cell->latmax[2] - M_PI_2;
/* And now the same, but with the longitudes... */
sctr_cell->lonmin[0] = 0.123456789;
sctr_cell->lonmax[0] = sctr_cell->lonmin[0] + 1.3235;
sctr_cell->lonmin[1] = 4.3445234;
sctr_cell->lonmax[1] = sctr_cell->lonmin[1] + 0.01;
sctr_cell->lonmin[2] = 0.0;
sctr_cell->lonmax[2] = sctr_cell->lonmin[2] + 2.0 * M_PI;
/* Finally, we define spherical radii of the scattered cells. Importantly,
* the spherical radius is *constants* over the cells (but may vary from
* cell to cell). */
sctr_cell->r[0] = shcs->r;
sctr_cell->r[1] = shcs->r + 1000.0;
sctr_cell->r[2] = shcs->r + 2000.0;
/* --------------------------------------------------------------------- */
/* Initialize an array to store the synthesized signal */
f = (double *)malloc(sctr_cell->ncell * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_cell(sctr_cell, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Print the synthesized values */
for (size_t i = 0; i < sctr_cell->ncell; i++)
printf("%0.16e\n", f[i]);
charm_crd_cell_free(sctr_cell);
free(f);
printf("===========================\n\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Solid spherical harmonic synthesis at a grid "
"of cells\n");
printf("===========================\n");
/* Initialize a "charm_cell" structure to hold a grid of cells with "15"
* cells in the latitudinal direction and "30" cells in the longitudinal
* direction. */
charm_cell *grd_cell = charm_crd_cell_calloc(CHARM_CRD_CELL_GRID, 15, 30);
/* Define some grid cells */
for (size_t i = 0; i < grd_cell->nlat; i++)
{
/* The maximum cell latitudes */
grd_cell->latmax[i] = M_PI_2 -
((double)i / (double)grd_cell->nlat) * M_PI;
/* The minimum cell latitudes */
grd_cell->latmin[i] = M_PI_2 -
((double)(i + 1) / (double)grd_cell->nlat) *
M_PI;
/* The spherical radii (may vary with the latitude, but we use
* a constant value here, because of the example that follows after
* this one). */
grd_cell->r[i] = shcs->r + 1000.0;
}
for (size_t j = 0; j < grd_cell->nlon; j++)
{
/* The minimum cell longitudes */
grd_cell->lonmin[j] = ((double)j / (double)grd_cell->nlon) *
(2.0 * M_PI);
/* The maximum cell longitudes */
grd_cell->lonmax[j] = ((double)(j + 1) / (double)grd_cell->nlon) *
(2.0 * M_PI);
}
/* Initialize an array to store the synthesized signal */
f = (double *)malloc(grd_cell->ncell * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_cell(grd_cell, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Print some of the synthesized values */
/* --------------------------------------------------------------------- */
i = 0;
j = 10;
printf("f(%zu, %zu) = %0.16e\n", i, j, f[i * grd_cell->nlon + j]);
i = 4;
j = 3;
printf("f(%zu, %zu) = %0.16e\n", i, j, f[i * grd_cell->nlon + j]);
/* --------------------------------------------------------------------- */
printf("===========================\n\n");
/* ===================================================================== */
/* ===================================================================== */
printf("\n");
printf("Spherical harmonic analysis of block-mean values in cells\n");
printf("===========================\n");
/* In this example, we use the "grd_cell" structure and the signal "f" from
* the previous example. */
/* Initialize a new structure of spherical harmonic coefficients for
* coefficients to be recovered from the harmonic analysis. Used are the
* same scalling constants as in "shcs". */
shcs2 = charm_shc_calloc(nmax, shcs->mu, shcs->r);
if (shcs2 == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.\n");
exit(CHARM_FAILURE);
}
/* Do the analysis using the method of the approximate quadrature */
charm_sha_cell(grd_cell, f, nmax, CHARM_SHA_CELL_AQ, shcs2, err);
charm_err_handler(err, 1);
/* Print some of the computed coefficients. Note that the harmonic
* analysis with block-mean values in cells in *not* exact, hence the
* coefficients will not be equal to the input ones. */
/* --------------------------------------------------------------------- */
n = 2;
unsigned long m = 0;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs2->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs2->s[m][n - m]);
n = 9;
m = 0;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs2->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs2->s[m][n - m]);
n = 9;
m = 4;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs2->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs2->s[m][n - m]);
n = 9;
m = 9;
printf("C(%3lu,%3lu) = %0.16e\n", n, m, shcs2->c[m][n - m]);
printf("S(%3lu,%3lu) = %0.16e\n", n, m, shcs2->s[m][n - m]);
/* --------------------------------------------------------------------- */
charm_crd_cell_free(grd_cell);
charm_shc_free(shcs2);
free(f);
printf("===========================\n\n");
/* ===================================================================== */
/* Free the heap memory */
/* ===================================================================== */
charm_err_free(err);
charm_shc_free(shcs);
printf("\nGreat, all done!\n");
exit(CHARM_SUCCESS);
/* ===================================================================== */
}
First- and second-order gradients (gravitational vector and tensor)#
This example shows how to compute the full gravitational vector and the full gravitational tensor in the Local north-oriented reference frame.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <charm/charm.h>
int main(void)
{
/* INPUTS */
/* ===================================================================== */
/* Path to the input file with spherical harmonic coefficients in the "gfc"
* format. */
char shcs_file[] = "../../data/input/EGM96-degree10.gfc";
/* ===================================================================== */
/* ===================================================================== */
/* Initialize a "charm_err" structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* First, let's get the maximum harmonic degree stored in "shcs_file".
* This is especially useful to read all coefficients in "shcs_file"
* without knowing its maximum harmonic degree a priori. */
unsigned long nmax = charm_shc_read_gfc(shcs_file, CHARM_SHC_NMAX_MODEL,
NULL, NULL, err);
charm_err_handler(err, 1);
/* Allocate memory for the spherical harmonic coefficients up to degree
* "nmax" */
charm_shc *shcs = charm_shc_malloc(nmax, 1.0, 1.0);
if (shcs == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_shc\" structure.\n");
exit(CHARM_FAILURE);
}
/* Now read all coefficients in "shcs_file" */
charm_shc_read_gfc(shcs_file, nmax, NULL, shcs, err);
charm_err_handler(err, 1);
/* Compute the Gauss--Legendre point grid for an "nmax" given by the
* maximum degree stored in "shcs" and for a radius equal to the reference
* sphere, to which the spherical harmonic coefficients are scaled to. We
* intentionally use here the Gauss--Legendre grid instead of the
* Driscoll--Healy grids in order to avoid the inaccuracies due to the
* singularities at the poles (see the documentation). */
charm_point *grd = charm_crd_point_gl(shcs->nmax, shcs->r);
if (grd == NULL)
{
fprintf(stderr, "Failed to compute the Gauss--Legendre grid.\n");
exit(CHARM_FAILURE);
}
/* Total number of gravitational vector elements */
size_t np = 3;
/* Allocate memory for the gravitational vector elements */
double **f = (double **)malloc(np * sizeof(double *));
if (f == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
for (size_t i = 0; i < np; i++)
{
f[i] = (double *)malloc(grd->npoint * sizeof(double));
if (f[i] == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
}
/* Compute the full first-order gradient (vector) */
charm_shs_point_grad1(grd, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Allocate memory for the magnitude of the gravitational vector */
double *g = (double *)malloc(grd->npoint * sizeof(double));
if (g == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
/* Compute the magnitude of the gravitational acceleration (no contribution
* due to the centrifugal force is considered here) */
for (size_t i = 0; i < grd->npoint; i++)
g[i] = sqrt(f[0][i] * f[0][i] + f[1][i] * f[1][i] + f[2][i] * f[2][i]);
/* Free the memory associated with the gravitational vector */
for (size_t i = 0; i < np; i++)
free(f[i]);
free(f);
/* Allocate memory for the gravitational tensor elements */
np = 6; /* There are six tensor elements */
f = (double **)malloc(np * sizeof(double *));
if (f == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
for (size_t i = 0; i < np; i++)
{
f[i] = (double *)malloc(grd->npoint * sizeof(double));
if (f[i] == NULL)
{
fprintf(stderr, "malloc failure\n");
exit(CHARM_FAILURE);
}
}
/* Now compute the full second-order gradient (tensor) */
charm_shs_point_grad2(grd, shcs, nmax, f, err);
charm_err_handler(err, 1);
/* Check whether "fxx + fyy + fzz" is zero within numerical errors */
double trace_error = fabs(f[0][0] + f[3][0] + f[5][0]); /* Trace at the
* first grid
* point of "grd"
* */
double tmp;
for (size_t i = 1; i < grd->npoint; i++)
{
tmp = fabs(f[0][i] + f[3][i] + f[5][i]);
if (tmp > trace_error)
trace_error = tmp;
}
printf("The largest error of the gravitational tensor trace is "
"%0.16e s^-2.\n", trace_error);
/* Free the heap memory */
charm_shc_free(shcs);
charm_crd_point_free(grd);
charm_err_free(err);
for (size_t i = 0; i < np; i++)
free(f[i]);
free(f);
free(g);
printf("Great, all done!\n");
/* ===================================================================== */
exit(CHARM_SUCCESS);
}
Fourier coefficients of Legendre functions#
#include <stdio.h>
#include <stdlib.h>
#include <charm/charm.h>
int main(void)
{
/* Maximum harmonic degree to compute the Fourier coefficients of Legendre
* functions */
unsigned long nmax = 500;
/* Initialize a structure to store the Fourier coefficients of Legendre
* functions */
charm_pnmj *pnmj = charm_leg_pnmj_calloc(nmax, CHARM_LEG_PMNJ);
if (pnmj == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_pnmj\" "
"structure.\n");
exit(CHARM_FAILURE);
}
/* Initialize an "charm_err" structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* Compute the Fourier coefficients */
charm_leg_pnmj_coeffs(pnmj, nmax, err);
charm_err_handler(err, 1);
/* Print some Fourier coefficients */
/* Harmonic degree */
unsigned long n = 123;
/* Harmonic order */
unsigned long m = 23;
/* Wave-number-related variable */
unsigned long j = 12;
/* Wave-number (computed from "j"; for details, see the "charm_leg"
* module)*/
unsigned long k = charm_leg_pnmj_j2k(n, j);
printf("Fourier coefficient for degree %lu, order %lu and "
"wave-number %lu = %0.16e\n", n, m, k, pnmj->pnmj[m][n - m][j]);
n = 360; m = 358; j = 101; k = charm_leg_pnmj_j2k(n, j);
printf("Fourier coefficient for degree %lu, order %lu and "
"wave-number %lu = %0.16e\n", n, m, k, pnmj->pnmj[m][n - m][j]);
/* Free the heap memory */
charm_leg_pnmj_free(pnmj);
charm_err_free(err);
printf("\nGreat, all done!\n");
exit(CHARM_SUCCESS);
}
Integrals#
Let’s compute an integral of a product of two spherical harmonics over a restricted domain on the unit sphere. Then, we do the same, but with the product of Legendre functions only.
#include <stdio.h>
#include <stdlib.h>
#include <charm/charm.h>
int main(void)
{
/* Type of the first spherical harmonic function, its harmonic degree and
* order, respectively */
_Bool i1 = 0; /* The "cos" spherical harmonic function */
unsigned long n1 = 287;
unsigned long m1 = 122;
/* Type of the second spherical harmonic function, its harmonic degree and
* order, respectively */
_Bool i2 = 1; /* The "sin" spherical harmonic function */
unsigned long n2 = 34;
unsigned long m2 = 9;
/* Minimum and maximum co-latitudes of the integration domain */
double cltmin = 0.1;
double cltmax = 0.9;
/* Minimum and maximum longitudes of the integration domain */
double lonmin = 0.5;
double lonmax = 0.6;
/* Initialize the Fourier coefficients of Legendre functions */
unsigned long nmax = CHARM_MAX(n1, n2);
charm_pnmj *pnmj = charm_leg_pnmj_calloc(nmax, CHARM_LEG_PMNJ);
if (pnmj == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_pnmj\" "
"structure.\n");
exit(CHARM_FAILURE);
}
/* Initialize an "charm_err" structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize the \"charm_err\" structure.\n");
exit(CHARM_FAILURE);
}
/* Compute the Fourier coefficients */
charm_leg_pnmj_coeffs(pnmj, nmax, err);
charm_err_handler(err, 1);
/* Compute the integral of a product of two spherical harmonics */
double iy = charm_integ_yi1n1m1yi2n2m2(cltmin, cltmax, lonmin, lonmax,
i1, n1, m1,
i2, n2, m2, pnmj, err);
charm_err_handler(err, 1);
/* Print the value of the integral */
printf("The integral of the product of two spherical harmonics "
"i1 = %d, n1 = %lu, m1 = %lu, i2 = %d, n2 = %lu, m2 = %lu "
"is %0.16e\n", i1, n1, m1, i2, n2, m2, iy);
/* Now compute the integral of a product of two Legendre functions over a
* restricted domain */
double ip = charm_integ_pn1m1pn2m2(cltmin, cltmax, n1, m1, n2, m2, pnmj,
err);
charm_err_handler(err, 1);
printf("The integral of the product of two Legendre functions "
"n1 = %lu, m1 = %lu, n2 = %lu, m2 = %lu "
"is %0.16e\n", n1, m1, n2, m2, ip);
/* Free the heap memory */
charm_err_free(err);
charm_leg_pnmj_free(pnmj);
exit(CHARM_SUCCESS);
}
Distributed computing with MPI#
This example shows how to use CHarm with distributed-memory systems using MPI.
Compile the source code with mpicc or similar, e.g., by:
mpicc mpi.c -l:libcharm.a -lfftw3 -lm
Launch the program using three MPI processes (required by the program), e.g.,:
mpiexec -n 3 ./a.out
As a side note, you can run the program using an actual distributed-memory system but also using an ordinary shared-memory system having at least 3 CPUs.
#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <charm/charm.h>
int main(int argc, char *argv[])
{
/* Initialize MPI execution environment */
/* --------------------------------------------------------------------- */
MPI_Init(&argc, &argv);
/* Note that the previous "MPI_Init" call does not allow to combined MPI
* and OpenMP parallelization techniques. For best performance, however,
* it is strongly advised to combine indeed MPI and OpenMP.
*
* Briefly, the recommended strategy is to use one MPI process per
* shared-memory system. This will minimize the data transfer between the
* computing nodes, which, in case of high-degree spherical harmonics
* transforms, is usually the bottleneck. Then, within each shared-memory
* system, all OpenMP threads can be used.
*
* For best performance, we recommend to use the MPI's
* "MPI_THREAD_FUNNELED" parallelization mode (only the thread that called
* "MPI_Init_thread" will make MPI calls). Remember that to combine MPI
* and OpenMP, the MPI execution environment has to be initialized like
* this:
*
* int provided;
* MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &provided);
* if (provided != MPI_THREAD_FUNNELED)
* {
* fprintf(stderr, "MPI didn't provide MPI_THREAD_FUNNELED\n");
* exit(CHARM_FAILURE);
* }
*
* */
/* --------------------------------------------------------------------- */
/* Identify who I am and the number of MPI processes */
/* --------------------------------------------------------------------- */
/* MPI communicator */
MPI_Comm comm = MPI_COMM_WORLD;
int rank; /* This is who I am */
int size; /* This is the total number of MPI processes */
MPI_Comm_rank(comm, &rank);
MPI_Comm_size(comm, &size);
/* We assume this code is launched at 3 MPI processes. */
if (size != 3)
{
fprintf(stderr, "This program must be launched at 3 MPI processes.\n");
exit(CHARM_FAILURE);
}
/* --------------------------------------------------------------------- */
/* Initialize an error structure */
/* --------------------------------------------------------------------- */
/* Note that we do not call "charm_err_init" with distributed computing */
charm_err *err = charm_mpi_err_init(comm);
/* If the previous function call failed at any process, "err" will be
* "NULL" at all processes in "comm". */
if (err == NULL)
{
fprintf(stderr, "Failed to create the \"charm_err\" structure");
exit(CHARM_FAILURE);
}
/* --------------------------------------------------------------------- */
/* Make up some fake spherical harmonic coefficients up to degree "10" */
/* --------------------------------------------------------------------- */
/* Maximum harmonic degree of spherical harmonic coefficients */
unsigned long nmax = 10;
/* Each MPI process may hold zero, one or more chunks of spherical harmonic
* coefficients. A chunk of spherical harmonic coefficients denotes all
* spherical harmonic coefficients of orders "m1, m1 + 1, ..., m2".
*
* In this example, the MPI processes with ranks "0" and "1" will hold two
* chunks while only a single chunk will be assigned to rank "2". This is
* just to demonstrate that the number of chunks may vary across MPI
* processes.
*
* Not shown here, but it is absolutely OK for some processes to store zero
* chunks, that is, no spherical harmonic coefficients at all. For
* instance, if you have "100" MPI processes and you want to create the
* "charm_shc" structure up to degree "0", you assign one chunk starting
* and ending at order "0" to one process and zero chunks to all the
* remaining "99" processes. */
size_t local_nchunk;
if (rank == 2)
local_nchunk = 1;
else
local_nchunk = 2;
/* To split spherical harmonic coefficients among MPI processes (that is,
* to define chunks), we need an array of "2 * local_nchunk" elements. For
* ranks "0" and "1", we thus need four elements but only two elements for
* rank "2". For simplicity, we allocate here, however, an array of four
* elements for all processes which is fine.
*
* In the example below, we hard-coded the chunks for demonstration
* purposes. However, this is likely not the strategy you will use in
* practice, especially if the number of MPI processes is large. Instead,
* you'll probably fill "local_order" based on the "local_nchunk", "size"
* and "rank" variables or so. */
unsigned long local_order[4];
if (rank == 0)
{
/* The first chunk will consist of all spherical harmonic coefficients
* of order "0" */
local_order[0] = 0;
local_order[1] = 0;
/* The second chunk will hold all coefficients of orders from "1" to
* "2" */
local_order[2] = 1;
local_order[3] = 2;
}
else if (rank == 1)
{
/* The first chunk at another MPI process will consist of all
* coefficients of orders from "3" to "5" */
local_order[0] = 3;
local_order[1] = 5;
/* The second chunk consists of all coefficients of order "6" */
local_order[2] = 6;
local_order[3] = 6;
}
else if (rank == 2)
{
/* This process has one chunk only, starting at order "7" and ending at
* "10" */
local_order[0] = 7;
local_order[1] = 10;
}
/* Note that CHarm offers rather flexible distribution of spherical
* harmonic coefficients among MPI processes. This has many benefits. If
* one of your computing nodes has lower/higher RAM capacity than the other
* nodes, you can assign lower/higher number of spherical harmonic
* coefficients to this process. Or you can distribute all coefficients
* evenly among all processes to ensure the memory consumption on all nodes
* will be even and so on and so forth. */
/* Now allocate the memory for a distributed "charm_shc" structure. The
* structure will be distributed among all MPI processed in "comm" based on
* our "local_nchunk" and "local_order". */
charm_shc *shcs = charm_mpi_shc_malloc(nmax, 1.0, 1.0, local_nchunk,
local_order, comm, err);
/* At this point, we have to check that the previous function call was
* successful. We can do this either be verifying that "shcs" is not
* "NULL" and/or we can use the CHarm's error handler.
*
* Many things could go wrong when calling "charm_mpi_shc_malloc" and other
* routines from the "mpi" module. For instance, there might be not enough
* memory available at some computing node(s), the partitioning of
* spherical harmonic coefficients to chunks might be wrong (e.g., by
* leaving some gaps by not specifying some spherical harmonic order(s) in
* "local_order", or by introducing overlaps between some chunks), the MPI
* communicator "comm" might be invalid and so on. You should therefore
* always call an error handler. It will try to describe the cause of the
* error, hopefully allowing you to easily fix it. */
charm_err_handler(err, 1);
/* Now fill "shcs" with some fake spherical harmonic coefficients.
* Usually, you will read these coefficients from files on your hard drive.
* But for simplicity, we use some artificial coefficients in this example.
* */
for (size_t j = 0; j < shcs->local_nchunk; j++)
{
for (unsigned long m = shcs->local_order[2 * j];
m <= shcs->local_order[2 * j + 1]; m++)
{
for (unsigned long n = m; n <= shcs->nmax; n++)
{
shcs->c[m][n - m] = (double)(m + n);
shcs->s[m][n - m] = (double)(m * n);
}
}
}
/* Important: Never access coefficients that are not owned by the process
* that is trying to access the coefficients.
*
* For instance, all coefficients of order "8" are now owned by (stored at)
* MPI process with rank "2". This
*
* "shcs->c[8][8 - 8]"
*
* will thus fail badly at processes with ranks other than "2". */
/* --------------------------------------------------------------------- */
/* Now let's prepare a Gauss--Legendre grid that will be distributed among
* the MPI processes */
/* --------------------------------------------------------------------- */
/* Get the number of latitudes and longitudes of the Gauss--Legendre grid
* that correspond to maximum harmonic degree "nmax" */
size_t nlat, nlon;
charm_crd_point_gl_shape(nmax, &nlat, &nlon);
/* For our "nmax = 10", "nlat" will be "11". We will thus have the
* following latitudes in the Gauss--Legendre grid: "lat0, lat1, lat2, ...,
* lat10". */
/* Similarly as we specified chunks of spherical harmonic coefficients, now
* we have to split the Gauss--Legendre grid into latitudinal chunks.
* Naturally, the sum of "local_nlat" across all processes much match
* "nlat". This is automatically checked when calling
* "charm_mpi_crd_point_gl" below.
*
* Again, assigning no latitudinal chunks to some MPI process is fairly OK.
* Below, this would translate into assigning "0" to "local_nlat" at one or
* more processes. This would of course necessitate to increase
* "local_nlat" at other processes. */
size_t local_nlat, local_0_start;
if (rank == 0)
{
/* This MPI process will hold 2 latitudes, "lat1" and "lat9". */
local_nlat = 2;
local_0_start = 1;
}
else if (rank == 1)
{
/* This process also holds two latitudes, "lat0" and "lat10". We
* intentionally put "lat0" to the process of rank "1" to demonstrate
* that CHarm does not care about which chunk is stored by which MPI
* process. */
local_nlat = 2;
local_0_start = 0;
}
else if (rank == 2)
{
/* This process holds seven latitudes, "lat2, lat3, ..., lat8" */
local_nlat = 7;
local_0_start = 2;
}
/* Note that "lat0 = -lat10", "lat1 = -lat9", ..., that is Gauss--Legendre
* grids are symmetric with respect to the equator. In case of symmetric
* grids in general, each process must hold the latitudes at both
* hemisphere. The only exception is that MPI process, which holds the
* equator (if there is any). In that case, "local_nlat" must be odd.
* This happens with Gauss--Legendre grids and even "nmax". See the
* documentation to "charm_mpi_crd_point_gl". */
charm_point *glg = charm_mpi_crd_point_gl(nmax, shcs->r, local_nlat,
local_0_start, comm, err);
charm_err_handler(err, 1);
/* --------------------------------------------------------------------- */
/* Let's print the latitudes, spherical radii and integration weights in
* "glg". The longitudes are not printed as they are the same at all
* processes. The data will be printed in the rank order. */
/* --------------------------------------------------------------------- */
for (int r = 0; r < size; r++)
{
if (rank == r)
{
for (size_t i = 0; i < glg->local_nlat; i++)
printf("rank %d: glg->lat[%zu] = %23.16e, "
"glg->r[%zu] = %23.16e, glg->w[%zu] = %23.16e\n",
rank, i, glg->lat[i], i, glg->r[i], i, glg->w[i]);
}
MPI_Barrier(comm);
}
/* --------------------------------------------------------------------- */
/* At this point, we have spherical harmonic coefficients in "shcs" as well
* as computing points in "glg" distributed across MPI processes, so we can
* do some distributed spherical harmonic synthesis and analysis. */
/* Spherical harmonic synthesis */
/* --------------------------------------------------------------------- */
/* At this point, each process holds "glg->local_npoint" points. The sum
* of "glg->local_npoint" across all MPI processes thus yields the total
* number of grid points "glg->npoint". Let's now allocate memory for the
* local signal "f" that is to be synthesized by each MPI process */
double *f = (double *)malloc(glg->local_npoint * sizeof(double));
if (f == NULL)
{
fprintf(stderr, "malloc failure.\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis */
charm_shs_point(glg, shcs, shcs->nmax, f, err);
charm_err_handler(err, 1);
/* Note that "f" is now distributed among MPI processes. The distribution
* is given by "glg". */
/* --------------------------------------------------------------------- */
/* Spherical harmonic analysis */
/* --------------------------------------------------------------------- */
/* Now we want to compute by spherical harmonic analysis coefficients from
* "f", which is distributed. First, we need to allocate a new distributed
* "charm_shc" structure. For simplicity, the distribution will be the
* same as with "shcs" but this is not required. */
charm_shc *shcs_out = charm_mpi_shc_malloc(nmax, 1.0, 1.0, local_nchunk,
local_order, comm, err);
charm_err_handler(err, 1);
/* Now the analysis. Note that the number of locally stored points in
* "glg" ("glg->local_npoint") must match the number of elements in "f" and
* that the elements of "f" must be ordered properly. */
charm_sha_point(glg, f, shcs_out->nmax, shcs_out, err);
charm_err_handler(err, 1);
/* "shcs" and "shcs_out" are the same, as can easily be verified. Their
* differences should not exceed the order of, say, "10**-13" or so. We
* print the outputs first in order according to the rank of the process
* and then in order given by "local_order". */
for (int r = 0; r < size; r++)
{
if (rank == r)
{
for (size_t j = 0; j < shcs->local_nchunk; j++)
{
for (unsigned long m = shcs->local_order[2 * j];
m <= shcs->local_order[2 * j + 1]; m++)
{
for (unsigned long n = m; n <= shcs->nmax; n++)
{
printf("rank %d: Cref_{%2lu, %2lu} = %0.16e, "
"Cnew_{%2lu, %2lu} = %0.16e, diff = %0.16e\n",
rank, n, m, shcs->c[m][n - m], n, m,
shcs_out->c[m][n - m],
shcs->c[m][n - m] - shcs_out->c[m][n - m]);
printf("rank %d: Sref_{%2lu, %2lu} = %0.16e, "
"Snew_{%2lu, %2lu} = %0.16e, diff = %0.16e\n",
rank, n, m, shcs->s[m][n - m], n, m,
shcs_out->s[m][n - m],
shcs->s[m][n - m] - shcs_out->s[m][n - m]);
}
}
}
}
MPI_Barrier(comm);
}
/* --------------------------------------------------------------------- */
/* IMPORTANT NOTE */
/* --------------------------------------------------------------------- */
/* CHarm has some global parameters that should be tuned to get the best
* performance with MPI. For instance, there is the
* "charm_glob_sha_block_lat_multiplier" variable, which has enormous
* impact on the performance of "charm_sha_point" in real-world high-degree
* applications.
*
* For demonstration purposes, let's change its value now and repeat the
* analysis. This has no effect on the results of course.
*
* See the documentation of the "glob" module for this and other associated
* variables. */
charm_glob_sha_block_lat_multiplier = 100;
charm_sha_point(glg, f, shcs_out->nmax, shcs_out, err);
charm_err_handler(err, 1);
/* --------------------------------------------------------------------- */
/* Free the memory */
/* --------------------------------------------------------------------- */
free(f);
charm_shc_free(shcs);
charm_shc_free(shcs_out);
charm_crd_point_free(glg);
charm_err_free(err);
/* --------------------------------------------------------------------- */
/* --------------------------------------------------------------------- */
printf("Great, all done by process %d of %d!\n", rank, size);
MPI_Finalize();
/* --------------------------------------------------------------------- */
exit(CHARM_SUCCESS);
/* ===================================================================== */
}
CHarm in single and quadruple precision#
A few simple rules need to be obeyed in single and quadruple precision.
Make sure you have compiled CHarm in single or quadruple precision (see Installing), whichever you need.
Replace every occurrence of the
charm_*prefix withcharmf_*for single precision andcharmq_*for quadruple precision. This applies to your own code, including the header files, and the compilation of your program.If the prefix is written in capital letters, do not alter it. The
CHARM_*symbols are common across all precisions of CHarm.In your code, replace every
doublevariable entering CHarm functions byfloatfor single precision or by__float128for quadruple precision.To compile your code in quadruple precision, link the
libquadmathlibrary, that is, add-lquadmathbefore-lm(see Compilation on Linux).When compiling your program, use the
fftwf*orfftwq*prefix to link FFTW in single or quadruple precision, respectively.
Example compilation in single precision#
Assuming your working directory is cookbook/c, an example compilation
with the static library might look like this:
gcc -fopenmp shcsf.c -l:libcharmf.a -lfftw3f -lfftw3f_omp -lm
Compilation with the shared library:
gcc -fopenmp -Wl,-rpath -Wl,/usr/local/lib shcsf.c \
-lcharmf -lfftw3f -lfftw3f_omp -lm
Example compilation in quadruple precision#
Static library:
gcc -fopenmp shcsq.c -l:libcharmq.a -lfftw3q -lfftw3q_omp \
-lquadmath -lm
Shared library:
gcc -fopenmp -Wl,-rpath -Wl,/usr/local/lib shcsq.c \
-lcharmq -lfftw3q -lfftw3q_omp -lquadmath -lm