petapm.h File Reference

#include <pfft.h>
#include "powerspectrum.h"

Include dependency graph for petapm.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	Region
struct	Layout
struct	PetaPMPriv
struct	PetaPM
struct	PetaPMParticleStruct
struct	PetaPMFunctions
struct	PetaPMGlobalFunctions

Typedefs
typedef struct Region	PetaPMRegion
typedef struct PetaPMPriv	PetaPMPriv
typedef struct PetaPM	PetaPM
typedef void(*	petapm_transfer_func) (PetaPM *pm, int64_t k2, int kpos[3], pfft_complex *value)
typedef void(*	petapm_readout_func) (PetaPM *pm, int i, double *mesh, double weight)
typedef PetaPMRegion *(*	petapm_prepare_func) (PetaPM *pm, PetaPMParticleStruct *pstruct, void *data, int *Nregions)
typedef void *(*	petapm_malloc_func) (char *name, size_t *size)
typedef void *(*	petapm_mfree_func) (void *ptr)

Functions
void	petapm_module_init (int Nthreads)
void	petapm_init (PetaPM *pm, double BoxSize, double Asmth, int Nmesh, double G, MPI_Comm comm)
void	petapm_destroy (PetaPM *pm)
void	petapm_region_init_strides (PetaPMRegion *region)
void	petapm_force (PetaPM *pm, petapm_prepare_func prepare, PetaPMGlobalFunctions *global_functions, PetaPMFunctions *functions, PetaPMParticleStruct *pstruct, void *userdata)
PetaPMRegion *	petapm_force_init (PetaPM *pm, petapm_prepare_func prepare, PetaPMParticleStruct *pstruct, int *Nregions, void *userdata)
pfft_complex *	petapm_force_r2c (PetaPM *pm, PetaPMGlobalFunctions *global_functions)
void	petapm_force_c2r (PetaPM *pm, pfft_complex *rho_k, PetaPMRegion *regions, const int Nregions, PetaPMFunctions *functions)
void	petapm_force_finish (PetaPM *pm)
PetaPMRegion *	petapm_get_fourier_region (PetaPM *pm)
PetaPMRegion *	petapm_get_real_region (PetaPM *pm)
int	petapm_mesh_to_k (PetaPM *pm, int i)
int *	petapm_get_thistask2d (PetaPM *pm)
int *	petapm_get_ntask2d (PetaPM *pm)
pfft_complex *	petapm_alloc_rhok (PetaPM *pm)

Typedef Documentation

PetaPM

typedef struct PetaPM PetaPM

petapm_malloc_func

typedef void*(* petapm_malloc_func) (char *name, size_t *size)

Definition at line 108 of file petapm.h.

petapm_mfree_func

typedef void*(* petapm_mfree_func) (void *ptr)

Definition at line 109 of file petapm.h.

petapm_prepare_func

typedef PetaPMRegion*(* petapm_prepare_func) (PetaPM *pm, PetaPMParticleStruct *pstruct, void *data, int *Nregions)

Definition at line 91 of file petapm.h.

petapm_readout_func

typedef void(* petapm_readout_func) (PetaPM *pm, int i, double *mesh, double weight)

Definition at line 90 of file petapm.h.

petapm_transfer_func

typedef void(* petapm_transfer_func) (PetaPM *pm, int64_t k2, int kpos[3], pfft_complex *value)

Definition at line 89 of file petapm.h.

PetaPMPriv

typedef struct PetaPMPriv PetaPMPriv

PetaPMRegion

typedef struct Region PetaPMRegion

Function Documentation

petapm_alloc_rhok()

pfft_complex* petapm_alloc_rhok

(

PetaPM *

pm

)

Definition at line 50 of file petapm.c.

{
    pfft_complex * rho_k = (pfft_complex * ) mymalloc("PMrho_k", pm->priv->fftsize * sizeof(double));
    memset(rho_k, 0, pm->priv->fftsize * sizeof(double));
    return rho_k;
}

References PetaPMPriv::fftsize, mymalloc, and PetaPM::priv.

Referenced by displacement_fields().

Here is the caller graph for this function:

petapm_destroy()

void petapm_destroy

(

PetaPM *

pm

)

Definition at line 215 of file petapm.c.

{
    pfft_destroy_plan(pm->priv->plan_forw);
    pfft_destroy_plan(pm->priv->plan_back);
    MPI_Comm_free(&pm->priv->comm_cart_2d);
    myfree(pm->Mesh2Task[0]);
}

References PetaPMPriv::comm_cart_2d, PetaPM::Mesh2Task, myfree, PetaPMPriv::plan_back, PetaPMPriv::plan_forw, and PetaPM::priv.

Referenced by do_force_test(), main(), and run_gravity_test().

Here is the caller graph for this function:

petapm_force()

void petapm_force	(	PetaPM *	pm,
		petapm_prepare_func	prepare,
		PetaPMGlobalFunctions *	global_functions,
		PetaPMFunctions *	functions,
		PetaPMParticleStruct *	pstruct,
		void *	userdata
	)

Definition at line 353 of file petapm.c.

                         {
    int Nregions;
    PetaPMRegion * regions = petapm_force_init(pm, prepare, pstruct, &Nregions, userdata);
    pfft_complex * rho_k = petapm_force_r2c(pm, global_functions);
    if(functions)
        petapm_force_c2r(pm, rho_k, regions, Nregions, functions);
    myfree(rho_k);
    if(CPS->RegionInd)
        myfree(CPS->RegionInd);
    myfree(regions);
    petapm_force_finish(pm);
}

References CPS, functions, global_functions, myfree, petapm_force_c2r(), petapm_force_finish(), petapm_force_init(), petapm_force_r2c(), and PetaPMParticleStruct::RegionInd.

Referenced by glass_force(), and gravpm_force().

Here is the call graph for this function:

Here is the caller graph for this function:

petapm_force_c2r()

void petapm_force_c2r	(	PetaPM *	pm,
		pfft_complex *	rho_k,
		PetaPMRegion *	regions,
		const int	Nregions,
		PetaPMFunctions *	functions
	)

Definition at line 317 of file petapm.c.

{
 
    PetaPMFunctions * f = functions;
    for (f = functions; f->name; f ++) {
        petapm_transfer_func transfer = f->transfer;
        petapm_readout_func readout = f->readout;
 
        pfft_complex * complx = (pfft_complex *) mymalloc("PMcomplex", pm->priv->fftsize * sizeof(double));
        /* apply the greens function turn rho_k into potential in fourier space */
        pm_apply_transfer_function(pm, rho_k, complx, transfer);
        walltime_measure("/PMgrav/calc");
 
        double * real = (double * ) mymalloc2("PMreal", pm->priv->fftsize * sizeof(double));
        pfft_execute_dft_c2r(pm->priv->plan_back, complx, real);
        walltime_measure("/PMgrav/c2r");
        myfree(complx);
        /* read out the potential: this will copy and free real.*/
        layout_build_and_exchange_cells_to_local(pm, &pm->priv->layout, pm->priv->meshbuf, real);
        walltime_measure("/PMgrav/comm");
 
        pm_iterate(pm, readout, regions, Nregions);
        walltime_measure("/PMgrav/readout");
    }
    walltime_measure("/PMgrav/Misc");
 
}

References PetaPMPriv::fftsize, functions, PetaPMPriv::layout, layout_build_and_exchange_cells_to_local(), PetaPMPriv::meshbuf, myfree, mymalloc, mymalloc2, PetaPMFunctions::name, PetaPMPriv::plan_back, pm_apply_transfer_function(), pm_iterate(), PetaPM::priv, PetaPMFunctions::readout, PetaPMFunctions::transfer, and walltime_measure.

Referenced by displacement_fields(), and petapm_force().

Here is the call graph for this function:

Here is the caller graph for this function:

petapm_force_finish()

void petapm_force_finish

(

PetaPM *

pm

)

Definition at line 348 of file petapm.c.

                                      {
    layout_finish(&pm->priv->layout);
    myfree(pm->priv->meshbuf);
}

References PetaPMPriv::layout, layout_finish(), PetaPMPriv::meshbuf, myfree, and PetaPM::priv.

Referenced by displacement_fields(), and petapm_force().

Here is the call graph for this function:

Here is the caller graph for this function:

petapm_force_init()

PetaPMRegion* petapm_force_init	(	PetaPM *	pm,
		petapm_prepare_func	prepare,
		PetaPMParticleStruct *	pstruct,
		int *	Nregions,
		void *	userdata
	)

Definition at line 251 of file petapm.c.

                         {
    CPS = pstruct;
 
    *Nregions = 0;
    PetaPMRegion * regions = prepare(pm, pstruct, userdata, Nregions);
    pm_init_regions(pm, regions, *Nregions);
 
    walltime_measure("/PMgrav/Misc");
    pm_iterate(pm, put_particle_to_mesh, regions, *Nregions);
    walltime_measure("/PMgrav/cic");
 
    layout_prepare(pm, &pm->priv->layout, pm->priv->meshbuf, regions, *Nregions, pm->comm);
 
    walltime_measure("/PMgrav/comm");
    return regions;
}

References PetaPM::comm, CPS, PetaPMPriv::layout, layout_prepare(), PetaPMPriv::meshbuf, pm_init_regions(), pm_iterate(), PetaPM::priv, put_particle_to_mesh(), and walltime_measure.

Referenced by displacement_fields(), and petapm_force().

Here is the call graph for this function:

Here is the caller graph for this function:

petapm_force_r2c()

pfft_complex* petapm_force_r2c	(	PetaPM *	pm,
		PetaPMGlobalFunctions *	global_functions
	)

Definition at line 273 of file petapm.c.

          {
    /* call pfft rho_k is CFT of rho */
 
    /* this is because
     *
     * CFT = DFT * dx **3
     * CFT[rho] = DFT [rho * dx **3] = DFT[CIC]
     * */
    double * real = (double * ) mymalloc2("PMreal", pm->priv->fftsize * sizeof(double));
    memset(real, 0, sizeof(double) * pm->priv->fftsize);
    layout_build_and_exchange_cells_to_pfft(pm, &pm->priv->layout, pm->priv->meshbuf, real);
    walltime_measure("/PMgrav/comm2");
 
#ifdef DEBUG
    verify_density_field(pm, real, pm->priv->meshbuf, pm->priv->meshbufsize);
    walltime_measure("/PMgrav/Misc");
#endif
 
    pfft_complex * complx = (pfft_complex *) mymalloc("PMcomplex", pm->priv->fftsize * sizeof(double));
    pfft_execute_dft_r2c(pm->priv->plan_forw, real, complx);
    myfree(real);
 
    pfft_complex * rho_k = (pfft_complex * ) mymalloc2("PMrho_k", pm->priv->fftsize * sizeof(double));
 
    /*Do any analysis that may be required before the transfer function is applied*/
    petapm_transfer_func global_readout = global_functions->global_readout;
    if(global_readout)
        pm_apply_transfer_function(pm, complx, rho_k, global_readout);
    if(global_functions->global_analysis)
        global_functions->global_analysis(pm);
    /*Apply the transfer function*/
    petapm_transfer_func global_transfer = global_functions->global_transfer;
    pm_apply_transfer_function(pm, complx, rho_k, global_transfer);
    walltime_measure("/PMgrav/r2c");
 
    report_memory_usage("PetaPM");
 
    myfree(complx);
    return rho_k;
}

References PetaPMPriv::fftsize, PetaPMGlobalFunctions::global_analysis, global_functions, PetaPMGlobalFunctions::global_readout, PetaPMGlobalFunctions::global_transfer, PetaPMPriv::layout, layout_build_and_exchange_cells_to_pfft(), PetaPMPriv::meshbuf, PetaPMPriv::meshbufsize, myfree, mymalloc, mymalloc2, PetaPMPriv::plan_forw, pm_apply_transfer_function(), PetaPM::priv, report_memory_usage, and walltime_measure.

Referenced by petapm_force().

Here is the call graph for this function:

Here is the caller graph for this function:

petapm_get_fourier_region()

PetaPMRegion* petapm_get_fourier_region

(

PetaPM *

pm

)

Definition at line 64 of file petapm.c.

                                                      {
    return &pm->fourier_space_region;
}

References PetaPM::fourier_space_region.

Referenced by displacement_fields().

Here is the caller graph for this function:

petapm_get_ntask2d()

int* petapm_get_ntask2d

(

PetaPM *

pm

)

Definition at line 77 of file petapm.c.

                                     {
    return pm->NTask2d;
}

References PetaPM::NTask2d.

Referenced by idgen_init().

Here is the caller graph for this function:

petapm_get_real_region()

PetaPMRegion* petapm_get_real_region

(

PetaPM *

pm

)

Definition at line 67 of file petapm.c.

                                                   {
    return &pm->real_space_region;
}

References PetaPM::real_space_region.

petapm_get_thistask2d()

int* petapm_get_thistask2d

(

PetaPM *

pm

)

Definition at line 74 of file petapm.c.

                                        {
    return pm->ThisTask2d;
}

References PetaPM::ThisTask2d.

Referenced by idgen_init().

Here is the caller graph for this function:

petapm_init()

void petapm_init	(	PetaPM *	pm,
		double	BoxSize,
		double	Asmth,
		int	Nmesh,
		double	G,
		MPI_Comm	comm
	)

Definition at line 94 of file petapm.c.

{
    /* define the global long / short range force cut */
    pm->BoxSize = BoxSize;
    pm->Asmth = Asmth;
    pm->Nmesh = Nmesh;
    pm->G = G;
    pm->CellSize = BoxSize / Nmesh;
    pm->comm = comm;
 
    ptrdiff_t n[3] = {Nmesh, Nmesh, Nmesh};
    ptrdiff_t np[2];
 
    int ThisTask;
    int NTask;
 
    pm->Mesh2Task[0] = (int *) mymalloc2("Mesh2Task", 2*sizeof(int) * Nmesh);
    pm->Mesh2Task[1] = pm->Mesh2Task[0] + Nmesh;
 
    MPI_Comm_rank(comm, &ThisTask);
    MPI_Comm_size(comm, &NTask);
 
    /* try to find a square 2d decomposition */
    int i;
    int k;
    for(i = sqrt(NTask) + 1; i >= 0; i --) {
        if(NTask % i == 0) break;
    }
    np[0] = i;
    np[1] = NTask / i;
 
    message(0, "Using 2D Task mesh %td x %td \n", np[0], np[1]);
    if( pfft_create_procmesh_2d(comm, np[0], np[1], &pm->priv->comm_cart_2d) ){
        endrun(0, "Error: This test file only works with %td processes.\n", np[0]*np[1]);
    }
 
    int periods_unused[2];
    MPI_Cart_get(pm->priv->comm_cart_2d, 2, pm->NTask2d, periods_unused, pm->ThisTask2d);
 
    if(pm->NTask2d[0] != np[0]) abort();
    if(pm->NTask2d[1] != np[1]) abort();
 
    pm->priv->fftsize = 2 * pfft_local_size_dft_r2c_3d(n, pm->priv->comm_cart_2d,
           PFFT_TRANSPOSED_OUT,
           pm->real_space_region.size, pm->real_space_region.offset,
           pm->fourier_space_region.size, pm->fourier_space_region.offset);
 
    /*
     * In fourier space, the transposed array is ordered in
     * are in (y, z, x). The strides and sizes returned
     * from local size is in (Nx, Ny, Nz), hence we roll them once
     * so that the strides will give correct linear indexing for
     * integer coordinates given in order of (y, z, x).
     * */
 
#define ROLL(a, N, j) { \
    typeof(a[0]) tmp[N]; \
    ptrdiff_t k; \
    for(k = 0; k < N; k ++) tmp[k] = a[k]; \
    for(k = 0; k < N; k ++) a[k] = tmp[(k + j)% N]; \
    }
 
    ROLL(pm->fourier_space_region.offset, 3, 1);
    ROLL(pm->fourier_space_region.size, 3, 1);
 
#undef ROLL
 
    /* calculate the strides */
    petapm_region_init_strides(&pm->real_space_region);
    petapm_region_init_strides(&pm->fourier_space_region);
 
    /* planning the fft; need temporary arrays */
 
    double * real = (double * ) mymalloc("PMreal", pm->priv->fftsize * sizeof(double));
    pfft_complex * rho_k = (pfft_complex * ) mymalloc("PMrho_k", pm->priv->fftsize * sizeof(double));
    pfft_complex * complx = (pfft_complex *) mymalloc("PMcomplex", pm->priv->fftsize * sizeof(double));
 
    pm->priv->plan_forw = pfft_plan_dft_r2c_3d(
        n, real, rho_k, pm->priv->comm_cart_2d, PFFT_FORWARD,
        PFFT_TRANSPOSED_OUT | PFFT_ESTIMATE | PFFT_TUNE | PFFT_DESTROY_INPUT);
    pm->priv->plan_back = pfft_plan_dft_c2r_3d(
        n, complx, real, pm->priv->comm_cart_2d, PFFT_BACKWARD,
        PFFT_TRANSPOSED_IN | PFFT_ESTIMATE | PFFT_TUNE | PFFT_DESTROY_INPUT);
 
    myfree(complx);
    myfree(rho_k);
    myfree(real);
 
    /* now lets fill up the mesh2task arrays */
 
#if 0
    message(1, "ThisTask = %d (%td %td %td) - (%td %td %td)\n", ThisTask,
            pm->real_space_region.offset[0],
            pm->real_space_region.offset[1],
            pm->real_space_region.offset[2],
            pm->real_space_region.size[0],
            pm->real_space_region.size[1],
            pm->real_space_region.size[2]);
#endif
 
    int * tmp = (int *) mymalloc("tmp", sizeof(int) * Nmesh);
    for(k = 0; k < 2; k ++) {
        for(i = 0; i < Nmesh; i ++) {
            tmp[i] = 0;
        }
        for(i = 0; i < pm->real_space_region.size[k]; i ++) {
            tmp[i + pm->real_space_region.offset[k]] = pm->ThisTask2d[k];
        }
        /* which column / row hosts this tile? */
        /* FIXME: this is very inefficient */
        MPI_Allreduce(tmp, pm->Mesh2Task[k], Nmesh, MPI_INT, MPI_MAX, comm);
        /*
        for(i = 0; i < Nmesh; i ++) {
            message(0, "Mesh2Task[%d][%d] == %d\n", k, i, Mesh2Task[k][i]);
        }
        */
    }
    myfree(tmp);
}

References PetaPM::Asmth, PetaPM::BoxSize, PetaPM::CellSize, PetaPM::comm, PetaPMPriv::comm_cart_2d, endrun(), PetaPMPriv::fftsize, PetaPM::fourier_space_region, PetaPM::G, G, PetaPM::Mesh2Task, message(), myfree, mymalloc, mymalloc2, PetaPM::Nmesh, NTask, PetaPM::NTask2d, Region::offset, petapm_region_init_strides(), PetaPMPriv::plan_back, PetaPMPriv::plan_forw, PetaPM::priv, PetaPM::real_space_region, ROLL, Region::size, ThisTask, and PetaPM::ThisTask2d.

Referenced by gravpm_init_periodic(), and main().

Here is the call graph for this function:

Here is the caller graph for this function:

petapm_mesh_to_k()

int petapm_mesh_to_k	(	PetaPM *	pm,
		int	i
	)

Definition at line 70 of file petapm.c.

                                         {
    /*Return the position of this point on the Fourier mesh*/
    return i<=pm->Nmesh/2 ? i : (i-pm->Nmesh);
}

References PetaPM::Nmesh.

Referenced by pm_apply_transfer_function().

Here is the caller graph for this function:

petapm_module_init()

void petapm_module_init

(

int

Nthreads

)

Definition at line 82 of file petapm.c.

{
    pfft_init();
 
    pfft_plan_with_nthreads(Nthreads);
 
    /* initialize the MPI Datatype of pencil */
    MPI_Type_contiguous(sizeof(struct Pencil), MPI_BYTE, &MPI_PENCIL);
    MPI_Type_commit(&MPI_PENCIL);
}

References MPI_PENCIL.

Referenced by begrun(), main(), and setup_tree().

Here is the caller graph for this function:

petapm_region_init_strides()

void petapm_region_init_strides

(

PetaPMRegion *

region

)

Definition at line 813 of file petapm.c.

                                                       {
    int k;
    size_t rt = 1;
    for(k = 2; k >= 0; k --) {
        region->strides[k] = rt;
        rt = region->size[k] * rt;
    }
    region->totalsize = rt;
    region->buffer = NULL;
}

References Region::buffer, Region::size, Region::strides, and Region::totalsize.

Referenced by _prepare(), convert_node_to_region(), makeregion(), and petapm_init().

Here is the caller graph for this function: