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);
}
Gravity forward modelling#
Given spherical harmonic coefficients of the Moon’s topography and of its density (constant, lateral and 3D density), this example computes the gravitational field implied by the lunar topographic masses.
The lunar shape is here defined by the MoonTopo2600p.shape model due to
Wieczorek, M. A., Gravity and topography of the terrestrial planets, Treatise on Geophysics, 10, 153-193, doi:10.1016/B978-0-444-53802-4.00169-X, 2015.
The density model is due to
Goossens, S., Sabaka, T.J., Wieczorek, M.A., Neumann, G.A., Mazarico, E., Lemoine, F., Nicholas, J.B., Smith, D.E., Zuber, M.T. (2020). High-resolution gravity field models from GRAIL data and implications for models of the density structure of the Moon’s crust, Journal of Geophysical Research: Planets, 125, e2019JE006086, doi:10.1029/2019JE006086.
The density model of Goossens et al. (2020) was transformed into the form required by our forward modelling method. While the transformed models that we use in this example are realistic, it is not recommended to use them in any real task. They were create for learning purposes only.
/* IMPORTANT: To compile this code, CHarm must be compiled with the MPFR
* support, that is, the "--enable-mpfr" installation flag must be used during
* the installation. If you didn't enable MPFR, manually delete (or comment
* out) all lines between the following marks:
*
* @@@
*
* lines to be deleted or commented out
*
* !!!
*
* */
#include <stdio.h>
#include <stdlib.h>
/* @@@ */
#include <mpfr.h>
/* !!! */
#include <charm/charm.h>
int main(void)
{
/* INPUTS */
/* ===================================================================== */
/* Shape of the gravitating body */
/* --------------------------------------------------------------------- */
/* Upper integration limit in spherical radius (topography) */
/* ..................................................................... */
/* Path to spherical harmonic coefficients defining the shape of the Moon
* */
char path_shcs_shape[] = "../../data/input/MoonTopo2600p_to10-tbl.txt";
/* Maximum harmonic degree of the lunar topography defined by coefficients
* in "path_shcs_shape" */
unsigned long shape_nmax = 10;
/* ..................................................................... */
/* Lower integration limit in spherical radius (sphere) */
/* ..................................................................... */
/* Radius of the reference sphere */
double rref = 1728000.0;
/* ..................................................................... */
/* --------------------------------------------------------------------- */
/* Density of the gravitating body */
/* --------------------------------------------------------------------- */
/* Path to spherical harmonic coefficients defining the zero- and
* first-order polynomial density coefficients */
char path_shcs_density0[] = "../../data/input/moon-rho0_to10.tbl";
char path_shcs_density1[] = "../../data/input/moon-rho1_to10.tbl";
char *path_shcs_density[2];
path_shcs_density[0] = path_shcs_density0;
path_shcs_density[1] = path_shcs_density1;
/* Maximum harmonic degrees of the polynomial density coefficients in
* "path_shcs_density". */
unsigned long density_nmax[2] = {10, 10};
/* Order of the density polynomial */
unsigned density_order = 1;
/* --------------------------------------------------------------------- */
/* Other inputs */
/* --------------------------------------------------------------------- */
/* Newton's gravitational constant ("kg^-1 * m^3 * s^-2") */
double gc = 6.67430e-11;
/* Mass of the Moon ("kg") */
double mass = 7.346e22;
/* Minimum and maximum topography powers. In theory, "pmax" should be
* infinitely large, so higher "pmax" means more complete forward modelling
* (but also longer computation times). */
unsigned pmin = 1;
unsigned pmax = 10;
/* Maximum harmonic degree of the output gravitational potential. Due to
* the non-linear relation between topography and its implied gravitational
* field, the potential series is spectrally unlimited even if "shape_nmax"
* is finite (the only exception is "shape_nmax = 0"). Therefore, this
* value should be as high as reasonably possible in order to mitigate
* forward modelling errors. */
unsigned long nmax_potential = 100;
/* --------------------------------------------------------------------- */
/* ===================================================================== */
/* ===================================================================== */
/* --------------------------------------------------------------------- */
/* Create CHarm's error structure */
charm_err *err = charm_err_init();
if (err == NULL)
{
fprintf(stderr, "Failed to initialize \"err\".\n");
exit(CHARM_FAILURE);
}
/* Allocate memory for spherical harmonic coefficients that will hold the
* shape coefficients */
charm_shc *shape_shcs = charm_shc_malloc(shape_nmax, 1.0, 1.0);
if (shape_shcs == NULL)
{
fprintf(stderr, "Failed to initialize \"shape_shcs\".\n");
exit(CHARM_FAILURE);
}
/* Read the shape coefficients */
charm_shc_read_tbl(path_shcs_shape, shape_nmax, shape_shcs, err);
charm_err_handler(err, 1);
/* Allocate memory for spherical harmonic coefficients that will hold the
* polynomial density coefficients and read the coefficients from input
* files */
charm_shc *density_shcs[2];
for (unsigned i = 0; i <= density_order; i++)
{
density_shcs[i] = charm_shc_malloc(density_nmax[i], 1.0, 1.0);
if (density_shcs[i] == NULL)
{
fprintf(stderr, "Failed to initialize \"density_shcs[i]\".\n");
exit(CHARM_FAILURE);
}
charm_shc_read_tbl(path_shcs_density[i], density_nmax[i],
density_shcs[i], err);
charm_err_handler(err, 1);
}
/* Allocate coefficients for the output gravitational potential */
charm_shc *potential_global_shcs = charm_shc_malloc(nmax_potential,
1.0, 1.0);
if (potential_global_shcs == NULL)
{
fprintf(stderr, "Failed to initialize \"potential_global_shcs\".\n");
exit(CHARM_FAILURE);
}
/* --------------------------------------------------------------------- */
/* Global gravity forward modelling using 3D density */
/* --------------------------------------------------------------------- */
charm_gfm_global_density_3d(shape_shcs, shape_nmax, rref,
density_shcs, density_nmax, density_order,
gc, mass,
pmin, pmax,
nmax_potential,
NULL, NULL, NULL,
potential_global_shcs,
err);
charm_err_handler(err, 1);
/* Now we can compute, say, the gravitational potential using
* "potential_global_shcs". Let's do this in the Gauss--Legendre grid that
* passes above all lunar masses. */
/* ..................................................................... */
/* Create the Gauss--Legendre grid */
charm_point *grd = charm_crd_point_gl(potential_global_shcs->nmax,
potential_global_shcs->r + 25000.0);
if (grd == NULL)
{
fprintf(stderr, "Failed to compute the Gauss--Legendre grid "
"\"grd\".\n");
exit(CHARM_FAILURE);
}
/* Allocate memory for the gravitational potential to be synthesized */
double *vgfm = (double *)malloc(grd->npoint * sizeof(double));
if (vgfm == NULL)
{
fprintf(stderr, "Failed to allocate memory for the gravitational "
"potential \"vgfm\".\n");
exit(CHARM_FAILURE);
}
/* Do the synthesis of the gravitational potential */
charm_shs_point(grd, potential_global_shcs, potential_global_shcs->nmax,
vgfm, err);
charm_err_handler(err, 1);
/* ..................................................................... */
/* --------------------------------------------------------------------- */
/* In CHarm, gravity forward modelling using lateral density is fairly
* similar to using 3D density, so we skip it. */
/* Global gravity forward modelling using constant density */
/* --------------------------------------------------------------------- */
/* Constant density of the lunar crust */
double density_const = 2550.0;
charm_gfm_global_density_const(shape_shcs, shape_nmax, rref,
density_const,
gc, mass,
pmin, pmax,
nmax_potential,
NULL, NULL, NULL,
potential_global_shcs,
err);
charm_err_handler(err, 1);
/* You could now similarly synthesize the gravitational potential, this
* time being implied by the constant density model. */
/* --------------------------------------------------------------------- */
/* @@@ */
/* Spatially restricted gravity forward modelling of near-zone masses using
* 3D density */
/* --------------------------------------------------------------------- */
/* Minimum and maximum order of the radial derivatives of the output
* quantity */
unsigned kmin = 0;
unsigned kmax = 3;
/* Radius of the sphere, on which evaluation points reside */
double r = rref + 25000.0;
/* Integration radius in radians */
double psi0 = 1.0;
/* Order of the potential derivative ("0" for potential, "1" for quantities
* related to gravitational vector elements, and "2" for quantities related
* to gravitational tensor elements) */
unsigned u = 0;
/* Order of the potential derivative with respect to the spherical distance
* */
unsigned v = 0;
/* We want to integrate all masses up to distance "psi0" from evalution
* points. To integrate masses beyond "psi0", set "zone" to
* "CHARM_GFM_FAR_ZONE". */
int zone = CHARM_GFM_NEAR_ZONE;
/* Number of bits to represent significands of floating points numbers used
* to evaluate truncation coefficients */
mpfr_prec_t nbits = 512;
charm_shc **potential_cap_shcs = malloc((kmax - kmin + 1) *
sizeof(charm_shc *));
if (potential_cap_shcs == NULL)
{
fprintf(stderr, "Failed to initialize \"potential_cap_shcs\".\n");
exit(CHARM_FAILURE);
}
for (unsigned k = kmin; k <= kmax; k++)
{
potential_cap_shcs[k - kmin] = charm_shc_malloc(nmax_potential,
1.0, 1.0);
if (potential_cap_shcs[k - kmin] == NULL)
{
fprintf(stderr, "Failed to initialize "
"\"potential_cap_shcs[k - kmin]\".\n");
exit(CHARM_FAILURE);
}
}
charm_gfm_cap_density_3d(shape_shcs, shape_nmax, rref,
density_shcs, density_nmax, density_order,
gc, mass,
pmin, pmax,
kmin, kmax,
r,
psi0,
u, v,
zone,
nbits,
nmax_potential,
NULL,
NULL,
NULL,
potential_cap_shcs,
err);
charm_err_handler(err, 1);
/* --------------------------------------------------------------------- */
/* Spatially restricted gravity forward modelling of near-zone masses using
* constant density */
/* --------------------------------------------------------------------- */
charm_gfm_cap_density_const(shape_shcs, shape_nmax, rref,
density_const,
gc, mass,
pmin, pmax,
kmin, kmax,
r,
psi0,
u, v,
zone,
nbits,
nmax_potential,
NULL,
NULL,
NULL,
potential_cap_shcs,
err);
charm_err_handler(err, 1);
/* --------------------------------------------------------------------- */
/* Let's compute some Molodensky's truncation coefficients, which are used
* internally whenever calling the "charm_gfm_cap_density_*" routines. */
/* --------------------------------------------------------------------- */
/* Allocate memory for trunction coefficients */
/* ..................................................................... */
/* Get the number of trunction coefficients */
size_t q_size = charm_gfm_cap_nq(nmax_potential, pmax, kmin, kmax,
density_order, err);
charm_err_handler(err, 1);
mpfr_t *q = (mpfr_t *)malloc(q_size * sizeof(mpfr_t));
if (q == NULL)
{
fprintf(stderr, CHARM_ERR_MALLOC_FAILURE);
exit(CHARM_FAILURE);
}
for (size_t j = 0; j < q_size; j++)
mpfr_init2(q[j], nbits);
/* ..................................................................... */
/* To compute the truncation coefficients, CHarm requires some floating
* point numbers to be provided as the "mpfr_t" data type, so let's do the
* conversion. We chose "rref", "r" and "psi0" such that they can be
* represented as "double"s exactly, so using "mpfr_set_d" should be safe.
* Otherwise, say, if "psi0 = 0.1234567", use "mpfr_set_str" instead. */
/* ..................................................................... */
mpfr_t mpfr_rref, mpfr_r, mpfr_psi0;
mpfr_inits2(nbits, mpfr_rref, mpfr_r, mpfr_psi0, (mpfr_ptr)NULL);
mpfr_set_d(mpfr_rref, rref, MPFR_RNDN);
mpfr_set_d(mpfr_r, r, MPFR_RNDN);
mpfr_set_d(mpfr_psi0, psi0, MPFR_RNDN);
/* ..................................................................... */
/* Compute the trucation coefficients */
/* ..................................................................... */
charm_gfm_cap_q(mpfr_rref, mpfr_r, mpfr_psi0,
nmax_potential, pmax, kmin, kmax, density_order,
CHARM_GFM_NEAR_ZONE, CHARM_GFM_Q00, nbits, q, err);
charm_err_handler(err, 1);
/* ..................................................................... */
/* Free the memory associated with the trunction coefficients */
/* ..................................................................... */
mpfr_clears(mpfr_rref, mpfr_r, mpfr_psi0, (mpfr_ptr)NULL);
for (size_t j = 0; j < q_size; j++)
mpfr_clear(q[j]);
free(q);
/* NOTE: At this point, the memory associated with "q" should be released.
* If you are linking against "glibc" (which is the case on most Linux
* distributions), the memory may not, however, be released completely.
* This is not because of some bug in MPFR or CHarm, but instead it seems
* to be the default behaviour of the GNU's "glibc" library. It seems that
* if you allocate a large amount of small memory chunks using malloc, as
* it happens in our code with
*
* for (size_t j = 0; j < q_size; j++)
* mpfr_init2(q[j], nbits);
*
* then "glibc" may not release the memory completely when we called
*
* for (size_t j = 0; j < q_size; j++)
* mpfr_clear(q[j]);
*
* In our example, "q_size" is fairly low though, so this should now be
* a problem in this case.
*
* If you are linking against "glibc", your "q_size" is large and you want
* to release all the memory, call now
*
* malloc_trim(0);
*
* Importantly, "malloc_trim" is a GNU extension, so you can call it only
* if you are linking agains "glibc". If you are using, say, "musl", you
* must not call "malloc_trim".
*
* Luckily, CHarm detects on compilation whether we are linking against
* "glibc". If so, "malloc_trim" is called before returning from functions
* dealing internally with truncation coefficients. Therefore, you do not
* have to worry that CHarm will leave some unreleased memory. But if you
* allocate memory for truncation coefficients outside CHarm, as we did in
* this example, it is your own responsibility to release or to not release
* the memory you allocated.
*
* For more details, look inside "src/mpfr/mpfr_flush_unreleased_memory.h".
*
* */
/* ..................................................................... */
/* --------------------------------------------------------------------- */
/* !!! */
/* Free the heap memory */
/* --------------------------------------------------------------------- */
charm_err_free(err);
charm_shc_free(shape_shcs);
for (unsigned i = 0; i <= density_order; i++)
charm_shc_free(density_shcs[i]);
charm_shc_free(potential_global_shcs);
for (unsigned k = kmin; k <= kmax; k++)
charm_shc_free(potential_cap_shcs[k - kmin]);
free(potential_cap_shcs);
charm_crd_point_free(grd);
free(vgfm);
/* --------------------------------------------------------------------- */
/* --------------------------------------------------------------------- */
printf("\nGreat, 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 indeed combine 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 must 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
* hemispheres. 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