petapm.c File Reference

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "types.h"
#include "petapm.h"
#include "utils.h"
#include "walltime.h"

Include dependency graph for petapm.c:

Go to the source code of this file.

Classes
struct	Pencil

Macros
#define	POS(i)   ((double*) (&((char*)CPS->Parts)[CPS->elsize * (i) + CPS->offset_pos]))
#define	MASS(i)   ((float*) (&((char*)CPS->Parts)[CPS->elsize * (i) + CPS->offset_mass]))
#define	INACTIVE(i)   (CPS->active && !CPS->active(i))
#define	ROLL(a, N, j)

Typedefs
typedef void(*	cell_iterator) (double *cell_value, double *comm_buffer)
typedef void(*	pm_iterator) (PetaPM *pm, int i, double *mesh, double weight)

Functions
static void	layout_prepare (PetaPM *pm, struct Layout *L, double *meshbuf, PetaPMRegion *regions, const int Nregions, MPI_Comm comm)
static void	layout_finish (struct Layout *L)
static void	layout_build_and_exchange_cells_to_pfft (PetaPM *pm, struct Layout *L, double *meshbuf, double *real)
static void	layout_build_and_exchange_cells_to_local (PetaPM *pm, struct Layout *L, double *meshbuf, double *real)
static void	layout_iterate_cells (PetaPM *pm, struct Layout *L, cell_iterator iter, double *real)
static int	pencil_cmp_target (const void *v1, const void *v2)
static int	pos_get_target (PetaPM *pm, const int pos[2])
static int64_t	reduce_int64 (int64_t input, MPI_Comm comm)
pfft_complex *	petapm_alloc_rhok (PetaPM *pm)
static void	pm_init_regions (PetaPM *pm, PetaPMRegion *regions, const int Nregions)
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)
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)
static void	pm_iterate (PetaPM *pm, pm_iterator iterator, PetaPMRegion *regions, const int Nregions)
static void	pm_apply_transfer_function (PetaPM *pm, pfft_complex *src, pfft_complex *dst, petapm_transfer_func H)
static void	put_particle_to_mesh (PetaPM *pm, int i, double *mesh, double weight)
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)
void	petapm_force (PetaPM *pm, petapm_prepare_func prepare, PetaPMGlobalFunctions *global_functions, PetaPMFunctions *functions, PetaPMParticleStruct *pstruct, void *userdata)
static void	layout_build_pencils (PetaPM *pm, struct Layout *L, double *meshbuf, PetaPMRegion *regions, const int Nregions)
static void	layout_exchange_pencils (struct Layout *L)
static void	to_pfft (double *cell, double *buf)
static void	to_region (double *cell, double *region)
static void	pm_iterate_one (PetaPM *pm, int i, pm_iterator iterator, PetaPMRegion *regions, const int Nregions)
void	petapm_region_init_strides (PetaPMRegion *region)

Variables
static MPI_Datatype	MPI_PENCIL
static PetaPMParticleStruct *	CPS

Macro Definition Documentation

INACTIVE

#define INACTIVE

(

i

)

(CPS->active && !CPS->active(i))

Definition at line 62 of file petapm.c.

MASS

#define MASS

(

i

)

((float*) (&((char*)CPS->Parts)[CPS->elsize * (i) + CPS->offset_mass]))

Definition at line 61 of file petapm.c.

POS

#define POS

(

i

)

((double*) (&((char*)CPS->Parts)[CPS->elsize * (i) + CPS->offset_pos]))

Definition at line 60 of file petapm.c.

ROLL

#define ROLL	(		a,
			N,
			j
	)

Value:

    { \
    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]; \
    }

Typedef Documentation

cell_iterator

typedef void(* cell_iterator) (double *cell_value, double *comm_buffer)

Definition at line 26 of file petapm.c.

pm_iterator

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

Definition at line 227 of file petapm.c.

Function Documentation

layout_build_and_exchange_cells_to_local()

static void layout_build_and_exchange_cells_to_local ( PetaPM *  pm, struct Layout *  L, double *  meshbuf, double *  real  )

static

Definition at line 639 of file petapm.c.

{
    L->BufRecv = (double *) mymalloc("PMBufRecv", L->NcImport * sizeof(double));
    int i;
    int offset;
 
    /*layout_iterate_cells transfers real to L->BufRecv*/
    layout_iterate_cells(pm, L, to_region, real);
 
    /*Real is done now: reuse the memory for BufSend*/
    myfree(real);
    /*Now allocate BufSend, which is confusingly used to receive data*/
    L->BufSend = (double *) mymalloc("PMBufSend", L->NcExport * sizeof(double));
 
    /* exchange cells */
    /* notice the order is reversed from to_pfft */
    MPI_Alltoallv(
            L->BufRecv, L->NcRecv, L->DcRecv, MPI_DOUBLE,
            L->BufSend, L->NcSend, L->DcSend, MPI_DOUBLE,
            L->comm);
 
    /* distribute BufSend to meshbuf */
    offset = 0;
    for(i = 0; i < L->NpExport; i ++) {
        struct Pencil * p = &L->PencilSend[i];
        memcpy(&meshbuf[p->meshbuf_first],
                L->BufSend + offset,
                sizeof(double) * p->len);
        offset += p->len;
    }
    myfree(L->BufSend);
    myfree(L->BufRecv);
}

References Layout::BufRecv, Layout::BufSend, Layout::comm, Layout::DcRecv, Layout::DcSend, layout_iterate_cells(), Pencil::len, Pencil::meshbuf_first, myfree, mymalloc, Layout::NcExport, Layout::NcImport, Layout::NcRecv, Layout::NcSend, Layout::NpExport, Pencil::offset, Layout::PencilSend, and to_region().

Referenced by petapm_force_c2r().

Here is the call graph for this function:

Here is the caller graph for this function:

layout_build_and_exchange_cells_to_pfft()

static void layout_build_and_exchange_cells_to_pfft ( PetaPM *  pm, struct Layout *  L, double *  meshbuf, double *  real  )

static

Definition at line 583 of file petapm.c.

{
    L->BufSend = (double *) mymalloc("PMBufSend", L->NcExport * sizeof(double));
    L->BufRecv = (double *) mymalloc("PMBufRecv", L->NcImport * sizeof(double));
 
    int i;
    int offset;
 
    /* collect all cells into the send buffer */
    offset = 0;
    for(i = 0; i < L->NpExport; i ++) {
        struct Pencil * p = &L->PencilSend[i];
        memcpy(L->BufSend + offset, &meshbuf[p->meshbuf_first],
                sizeof(double) * p->len);
        offset += p->len;
    }
 
    /* receive cells */
    MPI_Alltoallv(
            L->BufSend, L->NcSend, L->DcSend, MPI_DOUBLE,
            L->BufRecv, L->NcRecv, L->DcRecv, MPI_DOUBLE,
            L->comm);
 
#if 0
    double massExport = 0;
    for(i = 0; i < L->NcExport; i ++) {
        massExport += L->BufSend[i];
    }
 
    double massImport = 0;
    for(i = 0; i < L->NcImport; i ++) {
        massImport += L->BufRecv[i];
    }
    double totmassExport;
    double totmassImport;
    MPI_Allreduce(&massExport, &totmassExport, 1, MPI_DOUBLE, MPI_SUM, L->comm);
    MPI_Allreduce(&massImport, &totmassImport, 1, MPI_DOUBLE, MPI_SUM, L->comm);
    message(0, "totmassExport = %g totmassImport = %g\n", totmassExport, totmassImport);
#endif
 
    layout_iterate_cells(pm, L, to_pfft, real);
    myfree(L->BufRecv);
    myfree(L->BufSend);
}

References Layout::BufRecv, Layout::BufSend, Layout::comm, Layout::DcRecv, Layout::DcSend, layout_iterate_cells(), Pencil::len, Pencil::meshbuf_first, message(), myfree, mymalloc, Layout::NcExport, Layout::NcImport, Layout::NcRecv, Layout::NcSend, Layout::NpExport, Pencil::offset, Layout::PencilSend, and to_pfft().

Referenced by petapm_force_r2c().

Here is the call graph for this function:

Here is the caller graph for this function:

layout_build_pencils()

static void layout_build_pencils ( PetaPM *  pm, struct Layout *  L, double *  meshbuf, PetaPMRegion *  regions, const int  Nregions  )

static

Definition at line 479 of file petapm.c.

{
    /* now build pencils to be exported */
    int p0 = 0;
    int r;
    for (r = 0; r < Nregions; r++) {
        int ix;
#pragma omp parallel for private(ix)
        for(ix = 0; ix < regions[r].size[0]; ix++) {
            int iy;
            for(iy = 0; iy < regions[r].size[1]; iy++) {
                int poffset = ix * regions[r].size[1] + iy;
                struct Pencil * p = &L->PencilSend[p0 + poffset];
 
                p->offset[0] = ix + regions[r].offset[0];
                p->offset[1] = iy + regions[r].offset[1];
                p->offset[2] = regions[r].offset[2];
                p->len = regions[r].size[2];
                p->meshbuf_first = (regions[r].buffer - meshbuf) +
                    regions[r].strides[0] * ix +
                    regions[r].strides[1] * iy;
                /* now lets compress the pencil */
                while((p->len > 0) && (meshbuf[p->meshbuf_first + p->len - 1] == 0.0)) {
                    p->len --;
                }
                while((p->len > 0) && (meshbuf[p->meshbuf_first] == 0.0)) {
                    p->len --;
                    p->meshbuf_first++;
                    p->offset[2] ++;
                }
 
                p->task = pos_get_target(pm, p->offset);
            }
        }
        p0 += regions[r].size[0] * regions[r].size[1];
    }
 
}

References Region::buffer, Pencil::len, Pencil::meshbuf_first, Pencil::offset, Region::offset, Layout::PencilSend, pos_get_target(), Region::size, and Pencil::task.

Referenced by layout_prepare().

Here is the call graph for this function:

Here is the caller graph for this function:

layout_exchange_pencils()

static void layout_exchange_pencils ( struct Layout *  L )

static

Definition at line 522 of file petapm.c.

                                                       {
    int i;
    int offset;
    int NTask;
    MPI_Comm_size(L->comm, &NTask);
    /* build the first pointers to refer to the correct relative buffer locations */
    /* note that the buffer hasn't bee assembled yet */
    offset = 0;
    for(i = 0; i < NTask; i ++) {
        int j;
        struct Pencil * p = &L->PencilSend[offset];
        if(L->NpSend[i] == 0) continue;
        p->first = 0;
        for(j = 1; j < L->NpSend[i]; j++) {
            p[j].first = p[j - 1].first + p[j - 1].len;
        }
        offset += L->NpSend[i];
    }
 
    MPI_Alltoallv(
            L->PencilSend, L->NpSend, L->DpSend, MPI_PENCIL,
            L->PencilRecv, L->NpRecv, L->DpRecv, MPI_PENCIL,
            L->comm);
 
    /* set first to point to absolute position in the full import cell buffer */
    offset = 0;
    for(i = 0; i < NTask; i ++) {
        struct Pencil * p = &L->PencilRecv[offset];
        int j;
        for(j = 0; j < L->NpRecv[i]; j++) {
            p[j].first += L->DcRecv[i];
        }
        offset += L->NpRecv[i];
    }
 
    /* set first to point to absolute position in the full export cell buffer */
    offset = 0;
    for(i = 0; i < NTask; i ++) {
        struct Pencil * p = &L->PencilSend[offset];
        int j;
        for(j = 0; j < L->NpSend[i]; j++) {
            p[j].first += L->DcSend[i];
        }
        offset += L->NpSend[i];
    }
}

References Layout::comm, Layout::DcRecv, Layout::DcSend, Layout::DpRecv, Layout::DpSend, Pencil::first, Pencil::len, MPI_PENCIL, Layout::NpRecv, Layout::NpSend, NTask, Pencil::offset, Layout::PencilRecv, and Layout::PencilSend.

Referenced by layout_prepare().

Here is the caller graph for this function:

layout_finish()

static void layout_finish ( struct Layout *  L )

static

Definition at line 569 of file petapm.c.

                                             {
    myfree(L->PencilRecv);
    myfree(L->PencilSend);
    myfree(L->ibuffer);
}

References Layout::ibuffer, myfree, Layout::PencilRecv, and Layout::PencilSend.

Referenced by petapm_force_finish().

Here is the caller graph for this function:

layout_iterate_cells()

static void layout_iterate_cells ( PetaPM *  pm, struct Layout *  L, cell_iterator  iter, double *  real  )

static

Definition at line 682 of file petapm.c.

{
    int i;
#pragma omp parallel for
    for(i = 0; i < L->NpImport; i ++) {
        struct Pencil * p = &L->PencilRecv[i];
        int k;
        ptrdiff_t linear0 = 0;
        for(k = 0; k < 2; k ++) {
            int ix = p->offset[k];
            while(ix < 0) ix += pm->Nmesh;
            while(ix >= pm->Nmesh) ix -= pm->Nmesh;
            ix -= pm->real_space_region.offset[k];
            if(ix >= pm->real_space_region.size[k]) {
                /* serious problem assumption about pfft layout was wrong*/
                endrun(1, "bad pfft: original k: %d ix: %d, cur ix: %d, region: off %d size %d\n", k, p->offset[k], ix, pm->real_space_region.offset[k], pm->real_space_region.size[k]);
            }
            linear0 += ix * pm->real_space_region.strides[k];
        }
        int j;
        for(j = 0; j < p->len; j ++) {
            int iz = p->offset[2] + j;
            while(iz < 0) iz += pm->Nmesh;
            while(iz >= pm->Nmesh) iz -= pm->Nmesh;
            if(iz >= pm->real_space_region.size[2]) {
                /* serious problem assmpution about pfft layout was wrong*/
                endrun(1, "bad pfft: original iz: %d, cur iz: %d, region: off %d size %d\n", p->offset[2], iz, pm->real_space_region.offset[2], pm->real_space_region.size[2]);
            }
            ptrdiff_t linear = iz * pm->real_space_region.strides[2] + linear0;
            /*
             * operate on the pencil, either modifying real or BufRecv
             * */
            iter(&real[linear], &L->BufRecv[p->first + j]);
        }
    }
}

References Layout::BufRecv, endrun(), Pencil::first, Pencil::len, PetaPM::Nmesh, Layout::NpImport, Pencil::offset, Region::offset, Layout::PencilRecv, PetaPM::real_space_region, Region::size, and Region::strides.

Referenced by layout_build_and_exchange_cells_to_local(), and layout_build_and_exchange_cells_to_pfft().

Here is the call graph for this function:

Here is the caller graph for this function:

layout_prepare()

static void layout_prepare ( PetaPM *  pm, struct Layout *  L, double *  meshbuf, PetaPMRegion *  regions, const int  Nregions, MPI_Comm  comm  )

static

Definition at line 375 of file petapm.c.

{
    int r;
    int i;
    int NTask;
    L->comm = comm;
 
    MPI_Comm_size(L->comm, &NTask);
 
    L->ibuffer = (int *) mymalloc("PMlayout", sizeof(int) * NTask * 8);
 
    memset(L->ibuffer, 0, sizeof(int) * NTask * 8);
    L->NpSend = &L->ibuffer[NTask * 0];
    L->NpRecv = &L->ibuffer[NTask * 1];
    L->NcSend = &L->ibuffer[NTask * 2];
    L->NcRecv = &L->ibuffer[NTask * 3];
    L->DcSend = &L->ibuffer[NTask * 4];
    L->DcRecv = &L->ibuffer[NTask * 5];
    L->DpSend = &L->ibuffer[NTask * 6];
    L->DpRecv = &L->ibuffer[NTask * 7];
 
    L->NpExport = 0;
    L->NcExport = 0;
    L->NpImport = 0;
    L->NcImport = 0;
 
    int NpAlloc = 0;
    /* count pencils until buffer would run out */
    for (r = 0; r < Nregions; r ++) {
        NpAlloc += regions[r].size[0] * regions[r].size[1];
    }
 
    L->PencilSend = (struct Pencil *) mymalloc("PencilSend", NpAlloc * sizeof(struct Pencil));
 
    layout_build_pencils(pm, L, meshbuf, regions, Nregions);
 
    /* sort the pencils by the target rank for ease of next step */
    qsort_openmp(L->PencilSend, NpAlloc, sizeof(struct Pencil), pencil_cmp_target);
    /* zero length pixels are moved to the tail */
 
    /* now shrink NpExport*/
    L->NpExport = NpAlloc;
    while(L->NpExport > 0 && L->PencilSend[L->NpExport - 1].len == 0) {
        L->NpExport --;
    }
 
    /* count total number of cells to be exported */
    int NcExport = 0;
    for(i = 0; i < L->NpExport; i++) {
        int task = L->PencilSend[i].task;
        L->NcSend[task] += L->PencilSend[i].len;
        NcExport += L->PencilSend[i].len;
        L->NpSend[task] ++;
    }
    L->NcExport = NcExport;
 
    MPI_Alltoall(L->NpSend, 1, MPI_INT, L->NpRecv, 1, MPI_INT, L->comm);
    MPI_Alltoall(L->NcSend, 1, MPI_INT, L->NcRecv, 1, MPI_INT, L->comm);
 
    /* build the displacement array; why doesn't MPI build these automatically? */
    L->DpSend[0] = 0; L->DpRecv[0] = 0;
    L->DcSend[0] = 0; L->DcRecv[0] = 0;
    for(i = 1; i < NTask; i ++) {
        L->DpSend[i] = L->NpSend[i - 1] + L->DpSend[i - 1];
        L->DpRecv[i] = L->NpRecv[i - 1] + L->DpRecv[i - 1];
        L->DcSend[i] = L->NcSend[i - 1] + L->DcSend[i - 1];
        L->DcRecv[i] = L->NcRecv[i - 1] + L->DcRecv[i - 1];
    }
    L->NpImport = L->DpRecv[NTask -1] + L->NpRecv[NTask -1];
    L->NcImport = L->DcRecv[NTask -1] + L->NcRecv[NTask -1];
 
    /* some checks */
    if(L->DpSend[NTask - 1] + L->NpSend[NTask -1] != L->NpExport) {
        endrun(1, "NpExport = %d NpSend=%d DpSend=%d\n", L->NpExport, L->NpSend[NTask -1], L->DpSend[NTask - 1]);
    }
    if(L->DcSend[NTask - 1] + L->NcSend[NTask -1] != L->NcExport) {
        endrun(1, "NcExport = %d NcSend=%d DcSend=%d\n", L->NcExport, L->NcSend[NTask -1], L->DcSend[NTask - 1]);
    }
    int64_t totNpAlloc = reduce_int64(NpAlloc, L->comm);
    int64_t totNpExport = reduce_int64(L->NpExport, L->comm);
    int64_t totNcExport = reduce_int64(L->NcExport, L->comm);
    int64_t totNpImport = reduce_int64(L->NpImport, L->comm);
    int64_t totNcImport = reduce_int64(L->NcImport, L->comm);
 
    if(totNpExport != totNpImport) {
        endrun(1, "totNpExport = %ld\n", totNpExport);
    }
    if(totNcExport != totNcImport) {
        endrun(1, "totNcExport = %ld\n", totNcExport);
    }
 
    /* exchange the pencils */
    message(0, "PetaPM:  %010ld/%010ld Pencils and %010ld Cells\n", totNpExport, totNpAlloc, totNcExport);
    L->PencilRecv = (struct Pencil *) mymalloc("PencilRecv", L->NpImport * sizeof(struct Pencil));
    memset(L->PencilRecv, 0xfc, L->NpImport * sizeof(struct Pencil));
    layout_exchange_pencils(L);
}

References Layout::comm, Layout::DcRecv, Layout::DcSend, Layout::DpRecv, Layout::DpSend, endrun(), Layout::ibuffer, layout_build_pencils(), layout_exchange_pencils(), Pencil::len, message(), mymalloc, Layout::NcExport, Layout::NcImport, Layout::NcRecv, Layout::NcSend, Layout::NpExport, Layout::NpImport, Layout::NpRecv, Layout::NpSend, NTask, pencil_cmp_target(), Layout::PencilRecv, Layout::PencilSend, qsort_openmp, reduce_int64(), Region::size, and Pencil::task.

Referenced by petapm_force_init().

Here is the call graph for this function:

Here is the caller graph for this function:

pencil_cmp_target()

static int pencil_cmp_target ( const void *  v1, const void *  v2  )

static

Definition at line 837 of file petapm.c.

                                                               {
    const struct Pencil * p1 = (const struct Pencil *) v1;
    const struct Pencil * p2 = (const struct Pencil *) v2;
    /* move zero length pixels to the end */
    if(p2->len == 0) return -1;
    if(p1->len == 0) return 1;
    int t1 = p1->task;
    int t2 = p2->task;
    return ((t2 < t1) - (t1 < t2)) * 2 +
        ((p2->meshbuf_first < p1->meshbuf_first) - (p1->meshbuf_first < p2->meshbuf_first));
}

References Pencil::len, Pencil::meshbuf_first, and Pencil::task.

Referenced by layout_prepare().

Here is the caller graph for this function:

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:

pm_apply_transfer_function()

static void pm_apply_transfer_function ( PetaPM *  pm, pfft_complex *  src, pfft_complex *  dst, petapm_transfer_func  H  )

static

Definition at line 883 of file petapm.c.

         {
    size_t ip = 0;
 
    PetaPMRegion * region = &pm->fourier_space_region;
 
#pragma omp parallel for
    for(ip = 0; ip < region->totalsize; ip ++) {
        ptrdiff_t tmp = ip;
        int pos[3];
        int kpos[3];
        int64_t k2 = 0.0;
        int k;
        for(k = 0; k < 3; k ++) {
            pos[k] = tmp / region->strides[k];
            tmp -= pos[k] * region->strides[k];
            /* lets get the abs pos on the grid*/
            pos[k] += region->offset[k];
            /* check */
            if(pos[k] >= pm->Nmesh) {
                endrun(1, "position didn't make sense\n");
            }
            kpos[k] = petapm_mesh_to_k(pm, pos[k]);
            /* Watch out the cast */
            k2 += ((int64_t)kpos[k]) * kpos[k];
        }
        /* swap 0 and 1 because fourier space was transposed */
        /* kpos is y, z, x */
        pos[0] = kpos[2];
        pos[1] = kpos[0];
        pos[2] = kpos[1];
        dst[ip][0] = src[ip][0];
        dst[ip][1] = src[ip][1];
        if(H) {
            H(pm, k2, pos, &dst[ip]);
        }
    }
 
}

References endrun(), PetaPM::fourier_space_region, PetaPM::Nmesh, Region::offset, petapm_mesh_to_k(), Region::strides, and Region::totalsize.

Referenced by petapm_force_c2r(), and petapm_force_r2c().

Here is the call graph for this function:

Here is the caller graph for this function:

pm_init_regions()

static void pm_init_regions ( PetaPM *  pm, PetaPMRegion *  regions, const int  Nregions  )

static

Definition at line 723 of file petapm.c.

{
    if(regions) {
        int i;
        size_t size = 0;
        for(i = 0 ; i < Nregions; i ++) {
            size += regions[i].totalsize;
        }
        pm->priv->meshbufsize = size;
        if ( size == 0 ) return;
        pm->priv->meshbuf = (double *) mymalloc("PMmesh", size * sizeof(double));
        /* this takes care of the padding */
        memset(pm->priv->meshbuf, 0, size * sizeof(double));
        size = 0;
        for(i = 0 ; i < Nregions; i ++) {
            regions[i].buffer = pm->priv->meshbuf + size;
            size += regions[i].totalsize;
        }
    }
}

References Region::buffer, PetaPMPriv::meshbuf, PetaPMPriv::meshbufsize, mymalloc, PetaPM::priv, and Region::totalsize.

Referenced by petapm_force_init().

Here is the caller graph for this function:

pm_iterate()

static void pm_iterate ( PetaPM *  pm, pm_iterator  iterator, PetaPMRegion *  regions, const int  Nregions  )

static

Definition at line 804 of file petapm.c.

                                                                                                      {
    int i;
#pragma omp parallel for
    for(i = 0; i < CPS->NumPart; i ++) {
        pm_iterate_one(pm, i, iterator, regions, Nregions);
    }
    MPIU_Barrier(pm->comm);
}

References PetaPM::comm, CPS, MPIU_Barrier, PetaPMParticleStruct::NumPart, and pm_iterate_one().

Referenced by petapm_force_c2r(), and petapm_force_init().

Here is the call graph for this function:

Here is the caller graph for this function:

pm_iterate_one()

static void pm_iterate_one ( PetaPM *  pm, int  i, pm_iterator  iterator, PetaPMRegion *  regions, const int  Nregions  )

static

Definition at line 746 of file petapm.c.

{
    int k;
    int iCell[3];  /* integer coordinate on the regional mesh */
    double Res[3]; /* residual*/
    double * Pos = POS(i);
    const int RegionInd = CPS->RegionInd ? CPS->RegionInd[i] : 0;
 
    /* Asserts that the swallowed particles are not considered (region -2).*/
    if(RegionInd < 0)
        return;
    /* This should never happen: it is pure paranoia and to avoid icc being crazy*/
    if(RegionInd >= Nregions)
        endrun(1, "Particle %d has region %d out of bounds %d\n", i, RegionInd, Nregions);
 
    PetaPMRegion * region = &regions[RegionInd];
    for(k = 0; k < 3; k++) {
        double tmp = Pos[k] / pm->CellSize;
        iCell[k] = floor(tmp);
        Res[k] = tmp - iCell[k];
        iCell[k] -= region->offset[k];
        /* seriously?! particles are supposed to be contained in cells */
        if(iCell[k] >= region->size[k] - 1 || iCell[k] < 0) {
            endrun(1, "particle out of cell better stop %d (k=%d) %g %g %g region: %td %td\n", iCell[k],k,
                Pos[0], Pos[1], Pos[2],
                region->offset[k], region->size[k]);
        }
    }
 
    int connection;
    for(connection = 0; connection < 8; connection++) {
        double weight = 1.0;
        size_t linear = 0;
        for(k = 0; k < 3; k++) {
            int offset = (connection >> k) & 1;
            int tmp = iCell[k] + offset;
            linear += tmp * region->strides[k];
            weight *= offset?
                /* offset == 1*/ (Res[k])    :
                /* offset == 0*/ (1 - Res[k]);
        }
        if(linear >= region->totalsize) {
            endrun(1, "particle linear index out of cell better stop\n");
        }
        iterator(pm, i, &region->buffer[linear], weight);
    }
}

References Region::buffer, PetaPM::CellSize, CPS, endrun(), Pencil::offset, Region::offset, POS, PetaPMParticleStruct::RegionInd, Region::size, Region::strides, and Region::totalsize.

Referenced by pm_iterate().

Here is the call graph for this function:

Here is the caller graph for this function:

pos_get_target()

static int pos_get_target ( PetaPM *  pm, const int  pos[2]  )

static

Definition at line 824 of file petapm.c.

                                                         {
    int k;
    int task2d[2];
    int rank;
    for(k = 0; k < 2; k ++) {
        int ix = pos[k];
        while(ix < 0) ix += pm->Nmesh;
        while(ix >= pm->Nmesh) ix -= pm->Nmesh;
        task2d[k] = pm->Mesh2Task[k][ix];
    }
    MPI_Cart_rank(pm->priv->comm_cart_2d, task2d, &rank);
    return rank;
}

References PetaPMPriv::comm_cart_2d, PetaPM::Mesh2Task, PetaPM::Nmesh, and PetaPM::priv.

Referenced by layout_build_pencils().

Here is the caller graph for this function:

put_particle_to_mesh()

static void put_particle_to_mesh ( PetaPM *  pm, int  i, double *  mesh, double  weight  )

static

Definition at line 929 of file petapm.c.

                                                                                   {
    double Mass = *MASS(i);
    if(INACTIVE(i))
        return;
#pragma omp atomic update
    mesh[0] += weight * Mass;
}

References INACTIVE, and MASS.

Referenced by petapm_force_init().

Here is the caller graph for this function:

reduce_int64()

static int64_t reduce_int64 ( int64_t  input, MPI_Comm  comm  )

static

Definition at line 936 of file petapm.c.

                                                          {
    int64_t result = 0;
    MPI_Allreduce(&input, &result, 1, MPI_INT64, MPI_SUM, comm);
    return result;
}

References MPI_INT64.

Referenced by layout_prepare().

Here is the caller graph for this function:

to_pfft()

static void to_pfft ( double *  cell, double *  buf  )

static

Definition at line 577 of file petapm.c.

                                                 {
#pragma omp atomic update
            cell[0] += buf[0];
}

Referenced by layout_build_and_exchange_cells_to_pfft().

Here is the caller graph for this function:

to_region()

static void to_region ( double *  cell, double *  region  )

static

Definition at line 634 of file petapm.c.

                                                      {
    *region = *cell;
}

Referenced by layout_build_and_exchange_cells_to_local().

Here is the caller graph for this function:

Variable Documentation

CPS

PetaPMParticleStruct* CPS

static

Definition at line 59 of file petapm.c.

Referenced by petapm_force(), petapm_force_init(), pm_iterate(), and pm_iterate_one().

MPI_PENCIL

MPI_Datatype MPI_PENCIL

static

Definition at line 46 of file petapm.c.

Referenced by layout_exchange_pencils(), and petapm_module_init().