Commit a5b63b7b authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

forgot to add these two files.


Former-commit-id: 3c295f0d3723b6a3b9bde4b958b72e7cdb9ccf26
parent 4739cede
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
/* Config parameters. */
#include "../config.h"
/* Some standard headers. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <limits.h>
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#endif
/* Local headers. */
#include "error.h"
#include "const.h"
#include "cycle.h"
#include "atomic.h"
#include "lock.h"
#include "space.h"
#include "part.h"
#include "multipole.h"
#include "cell.h"
/**
* @brief Merge two multipoles.
*
* @param ma The #multipole which will contain the merged result.
* @param mb The other #multipole.
*/
void multipole_merge ( struct multipole *ma , struct multipole *mb ) {
#if multipole_order == 1
/* Correct the position. */
float mma = ma->coeffs[0], mmb = mb->coeffs[0];
float w = 1.0f / ( mma + mmb );
for ( int k = 0 ; k < 3 ; k++ )
ma->x[k] = ( ma->x[k]*mma + mb->x[k]*mmb ) * w;
/* Add the particle to the moments. */
ma->coeffs[0] = mma + mmb;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Add a particle to the given multipole.
*
* @param m The #multipole.
* @param p The #gpart.
*/
void multipole_addpart ( struct multipole *m , struct gpart *p ) {
#if multipole_order == 1
/* Correct the position. */
float mm = m->coeffs[0], mp = p->mass;
float w = 1.0f / ( mm + mp );
for ( int k = 0 ; k < 3 ; k++ )
m->x[k] = ( m->x[k]*mm + p->x[k]*mp ) * w;
/* Add the particle to the moments. */
m->coeffs[0] = mm + mp;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Add a group of particles to the given multipole.
*
* @param m The #multipole.
* @param p The #gpart array.
* @param N Number of parts to add.
*/
void multipole_addparts ( struct multipole *m , struct gpart *p , int N ) {
#if multipole_order == 1
/* Get the combined mass and positions. */
double xp[3] = { 0.0 , 0.0 , 0.0 };
float mp = 0.0f, w;
for ( int k = 0 ; k < N ; k++ ) {
w = p[k].mass;
mp += w;
xp[0] += p[k].x[0] * w;
xp[1] += p[k].x[1] * w;
xp[2] += p[k].x[2] * w;
}
/* Correct the position. */
float mm = m->coeffs[0];
w = 1.0f / ( mm + mp );
for ( int k = 0 ; k < 3 ; k++ )
m->x[k] = ( m->x[k]*mm + xp[k] ) * w;
/* Add the particle to the moments. */
m->coeffs[0] = mm + mp;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Init a multipole from a set of particles.
*
* @param m The #multipole.
* @param parts The #gparts.
* @param N The number of particles.
*/
void multipole_init ( struct multipole *m , struct gpart *parts , int N ) {
#if multipole_order == 1
float mass = 0.0f, w;
double x[3] = { 0.0 , 0.0 , 0.0 };
int k;
/* Collect the particle data. */
for ( k = 0 ; k < N ; k++ ) {
w = parts[k].mass;
mass += w;
x[0] += parts[k].x[0] * w;
x[1] += parts[k].x[1] * w;
x[2] += parts[k].x[2] * w;
}
/* Store the data on the multipole. */
m->coeffs[0] = mass;
m->x[0] = x[0] / mass;
m->x[1] = x[1] / mass;
m->x[2] = x[2] / mass;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Reset the data of a #multipole.
*
* @param m The #multipole.
*/
void multipole_reset ( struct multipole *m ) {
/* Just bzero the struct. */
bzero( m , sizeof(struct multipole) );
}
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
/* Some constants. */
#define multipole_order 1
/* Multipole struct. */
struct multipole {
/* Multipole location. */
double x[3];
/* Acceleration on this multipole. */
float a[3];
/* Multipole coefficients. */
float coeffs[ multipole_order*multipole_order ];
};
/* Multipole function prototypes. */
static void multipole_iact_mm ( struct multipole *ma , struct multipole *mb , double *shift );
void multipole_merge ( struct multipole *ma , struct multipole *mb );
void multipole_addpart ( struct multipole *m , struct gpart *p );
void multipole_addparts ( struct multipole *m , struct gpart *p , int N );
void multipole_init ( struct multipole *m , struct gpart *parts , int N );
void multipole_reset ( struct multipole *m );
#include <math.h>
#include "kernel.h"
/**
* @brief Compute the pairwise interaction between two multipoles.
*
* @param ma The first #multipole.
* @param mb The second #multipole.
* @param shift The periodicity correction.
*/
__attribute__ ((always_inline)) INLINE static void multipole_iact_mm ( struct multipole *ma , struct multipole *mb , double *shift ) {
float dx[3], ir, r, r2 = 0.0f, acc;
int k;
/* Compute the multipole distance. */
for ( k = 0 ; k < 3 ; k++ ) {
dx[k] = ma->x[k] - mb->x[k] - shift[k];
r2 += dx[k]*dx[k];
}
/* Compute the normalized distance vector. */
ir = 1.0f / sqrtf( r2 );
r = r2 * ir;
/* Evaluate the gravity kernel. */
kernel_grav_eval( r , &acc );
/* Scale the acceleration. */
acc *= const_G * ir * ir * ir;
/* Compute the forces on both multipoles. */
#if multipole_order == 1
float mma = ma->coeffs[0], mmb = mb->coeffs[0];
for ( k = 0 ; k < 3 ; k++ ) {
ma->a[k] -= dx[k] * acc * mmb;
mb->a[k] += dx[k] * acc * mma;
}
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
/**
* @brief Compute the interaction of a multipole on a particle.
*
* @param m The #multipole.
* @param p The #gpart.
* @param shift The periodicity correction.
*/
__attribute__ ((always_inline)) INLINE static void multipole_iact_mp ( struct multipole *m , struct gpart *p , double *shift ) {
float dx[3], ir, r, r2 = 0.0f, acc;
int k;
/* Compute the multipole distance. */
for ( k = 0 ; k < 3 ; k++ ) {
dx[k] = m->x[k] - p->x[k] - shift[k];
r2 += dx[k]*dx[k];
}
/* Compute the normalized distance vector. */
ir = 1.0f / sqrtf( r2 );
r = r2 * ir;
/* Evaluate the gravity kernel. */
kernel_grav_eval( r , &acc );
/* Scale the acceleration. */
acc *= const_G * ir * ir * ir * m->coeffs[0];
/* Compute the forces on both multipoles. */
#if multipole_order == 1
for ( k = 0 ; k < 3 ; k++ )
p->a[k] += dx[k] * acc;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment