Skip to content
Snippets Groups Projects
Commit 3fd578a4 authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

major re-factoring of interactions and added separation between density and force computations. 

Former-commit-id: 30d7490ecc978a2a5a8961009ade54b61d5fa69b
parent 94fa733e
Branches
Tags
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
src/cell.h
+ 3
− 0
View file @ 3fd578a4
...@@ -64,6 +64,9 @@ struct cell { ...@@ -64,6 +64,9 @@ struct cell {
/* The tasks computing this cell's sorts. */ /* The tasks computing this cell's sorts. */
struct task *sorts[14]; struct task *sorts[14];
/* The ghost task to link density to interactions. */
struct task *ghost;
/* Number of tasks this cell is waiting for and whether it is in use. */ /* Number of tasks this cell is waiting for and whether it is in use. */
int wait; int wait;
......
src/gadgetsmp.h
+ 1
− 0
View file @ 3fd578a4
...@@ -28,5 +28,6 @@ ...@@ -28,5 +28,6 @@
#include "cell.h" #include "cell.h"
#include "space.h" #include "space.h"
#include "queue.h" #include "queue.h"
#include "runner_iact.h"
#include "runner.h" #include "runner.h"
src/part.h
+ 29
− 6
View file @ 3fd578a4
...@@ -26,11 +26,8 @@ ...@@ -26,11 +26,8 @@
/* Data of a single particle. */ /* Data of a single particle. */
struct part { struct part {
/* Particle position. */
double x[3];
/* Particle cutoff radius. */ /* Particle cutoff radius. */
float r; float h;
/* Particle time-step. */ /* Particle time-step. */
float dt; float dt;
...@@ -38,8 +35,34 @@ struct part { ...@@ -38,8 +35,34 @@ struct part {
/* Particle ID. */ /* Particle ID. */
int id; int id;
/* Number of pairwise interactions. */ /* Particle density. */
double count, count_dh; float rho;
/* Particle position. */
double x[3];
/* Particle velocity. */
float v[3];
/* Particle acceleration. */
float a[3];
/* Particle pressure. */
float P;
/* Particle mass. */
float m;
/* Particle internal energy. */
float u;
/* Change in particle energy over time. */
float u_dt;
/* Derivative of the density with respect to this particle's smoothing length. */
float rho_dh;
/* Particle number density. */
int icount; int icount;
} __attribute__((aligned (32))); } __attribute__((aligned (32)));
......
src/queue.c
+ 11
− 11
View file @ 3fd578a4
...@@ -185,11 +185,11 @@ struct task *queue_gettask ( struct queue *q , int blocking , int keep ) { ...@@ -185,11 +185,11 @@ struct task *queue_gettask ( struct queue *q , int blocking , int keep ) {
continue; continue;
/* Different criteria for different types. */ /* Different criteria for different types. */
if ( res->type == tid_self || (res->type == tid_sub && res->cj == NULL) ) { if ( res->type == task_type_self || (res->type == task_type_sub && res->cj == NULL) ) {
if ( res->ci->hold || cell_locktree( res->ci ) != 0 ) if ( res->ci->hold || cell_locktree( res->ci ) != 0 )
continue; continue;
} }
else if ( res->type == tid_pair || (res->type == tid_sub && res->cj != NULL) ) { else if ( res->type == task_type_pair || (res->type == task_type_sub && res->cj != NULL) ) {
if ( res->ci->hold || res->cj->hold || res->ci->wait || res->cj->wait ) if ( res->ci->hold || res->cj->hold || res->ci->wait || res->cj->wait )
continue; continue;
if ( cell_locktree( res->ci ) != 0 ) if ( cell_locktree( res->ci ) != 0 )
...@@ -306,7 +306,7 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int ...@@ -306,7 +306,7 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int
continue; continue;
/* Get the score for this task. */ /* Get the score for this task. */
if ( res->type == tid_self || res->type == tid_sort || ( res->type == tid_sub && res->cj == NULL ) ) if ( res->type == task_type_self || res->type == task_type_sort || ( res->type == task_type_sub && res->cj == NULL ) )
score = ( res->ci->owner == rid ); score = ( res->ci->owner == rid );
else else
score = ( res->ci->owner == rid ) + ( res->cj->owner == rid ); score = ( res->ci->owner == rid ) + ( res->cj->owner == rid );
...@@ -314,11 +314,11 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int ...@@ -314,11 +314,11 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int
continue; continue;
/* Different criteria for different types. */ /* Different criteria for different types. */
if ( res->type == tid_self || (res->type == tid_sub && res->cj == NULL) ) { if ( res->type == task_type_self || (res->type == task_type_sub && res->cj == NULL) ) {
if ( res->ci->hold || cell_locktree( res->ci ) != 0 ) if ( res->ci->hold || cell_locktree( res->ci ) != 0 )
continue; continue;
} }
else if ( res->type == tid_pair || (res->type == tid_sub && res->cj != NULL) ) { else if ( res->type == task_type_pair || (res->type == task_type_sub && res->cj != NULL) ) {
if ( res->ci->hold || res->cj->hold || res->ci->wait || res->cj->wait ) if ( res->ci->hold || res->cj->hold || res->ci->wait || res->cj->wait )
continue; continue;
if ( cell_locktree( res->ci ) != 0 ) if ( cell_locktree( res->ci ) != 0 )
...@@ -332,9 +332,9 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int ...@@ -332,9 +332,9 @@ struct task *queue_gettask_new ( struct queue *q , int rid , int blocking , int
/* If we owned a previous task, unlock it. */ /* If we owned a previous task, unlock it. */
if ( ind_best >= 0 ) { if ( ind_best >= 0 ) {
res = &qtasks[ qtid[ ind_best ] ]; res = &qtasks[ qtid[ ind_best ] ];
if ( res->type == tid_self || res->type == tid_pair || res->type == tid_sub ) if ( res->type == task_type_self || res->type == task_type_pair || res->type == task_type_sub )
cell_unlocktree( res->ci ); cell_unlocktree( res->ci );
if ( res->type == tid_pair || (res->type == tid_sub && res->cj != NULL) ) if ( res->type == task_type_pair || (res->type == task_type_sub && res->cj != NULL) )
cell_unlocktree( res->cj ); cell_unlocktree( res->cj );
} }
...@@ -438,19 +438,19 @@ void queue_sort ( struct queue *q ) { ...@@ -438,19 +438,19 @@ void queue_sort ( struct queue *q ) {
for ( k = 0 ; k < q->count ; k++ ) { for ( k = 0 ; k < q->count ; k++ ) {
t = &q->tasks[ q->tid[k] ]; t = &q->tasks[ q->tid[k] ];
switch ( t->type ) { switch ( t->type ) {
case tid_self: case task_type_self:
wait[k] = t->rank; wait[k] = t->rank;
weight[k] = 0; // t->ci->count * t->ci->count; weight[k] = 0; // t->ci->count * t->ci->count;
break; break;
case tid_pair: case task_type_pair:
wait[k] = t->rank; wait[k] = t->rank;
weight[k] = 0; // t->ci->count * t->cj->count; weight[k] = 0; // t->ci->count * t->cj->count;
break; break;
case tid_sub: case task_type_sub:
wait[k] = t->rank; wait[k] = t->rank;
weight[k] = 0; // (t->cj == NULL) ? t->ci->count * t->ci->count : t->ci->count * t->cj->count; weight[k] = 0; // (t->cj == NULL) ? t->ci->count * t->ci->count : t->ci->count * t->cj->count;
break; break;
case tid_sort: case task_type_sort:
wait[k] = t->rank; wait[k] = t->rank;
weight[k] = 0; // t->ci->count; weight[k] = 0; // t->ci->count;
break; break;
......
src/runner.c
+ 36
− 508
View file @ 3fd578a4
This diff is collapsed.
src/runner.h
+ 3
− 127
View file @ 3fd578a4
...@@ -109,131 +109,6 @@ long long int runner_hist_bins[ runner_hist_N ]; ...@@ -109,131 +109,6 @@ long long int runner_hist_bins[ runner_hist_N ];
#endif #endif
/**
* @brief Compute the 'interaction' between two particles.
*
* @param r2 the inter-particle radius squared.
* @param hi the screening distance of the ith particle.
* @param hj the screening distance of the jth particle.
* @param io Pointer to where to store the interaction of the ith particle.
* @param jo Pointer to where to store the interaction of the ith particle.
*/
__attribute__ ((always_inline)) INLINE void iact_nopart ( float r2 , float hi , float hj , float *force_i , float *force_j , int *count_i , int *count_j ) {
#define KERNEL_COEFF_1 2.546479089470f
#define KERNEL_COEFF_2 15.278874536822f
#define KERNEL_COEFF_3 45.836623610466f
#define KERNEL_COEFF_4 30.557749073644f
#define KERNEL_COEFF_5 5.092958178941f
#define KERNEL_COEFF_6 (-15.278874536822f)
#define NORM_COEFF 4.188790204786f
float r = sqrtf( r2 );
float ui, uj, wi, wj;
if ( r2 < hi*hi && !( force_i == NULL && count_i == NULL ) ) {
ui = r / hi;
if ( ui < 0.5 )
wi = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (ui - 1.0f) * ui * ui;
else
wi = KERNEL_COEFF_5 * (1.0f - ui) * (1.0f - ui) * (1.0 - ui);
if ( force_i != NULL )
*force_i += NORM_COEFF * wi;
if ( count_i != NULL )
*count_i += 1;
}
if ( r2 < hj*hj && !( force_j == NULL && count_j == NULL ) ) {
uj = r / hj;
if ( uj < 0.5 )
wj = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (uj - 1.0f) * uj * uj;
else
wj = KERNEL_COEFF_5 * (1.0f - uj) * (1.0f - uj) * (1.0 - uj);
if ( force_j != NULL )
*force_j += NORM_COEFF * wj;
if ( count_j != NULL )
*count_j += 1;
}
#ifdef HIST
if ( hi > hj )
runner_hist_hit( hi / hj );
else
runner_hist_hit( hj / hi );
#endif
}
__attribute__ ((always_inline)) INLINE void iact ( float r2 , float hi , float hj , struct part *pi , struct part *pj ) {
#define KERNEL_COEFF_1 2.546479089470f
#define KERNEL_COEFF_2 15.278874536822f
#define KERNEL_COEFF_3 45.836623610466f
#define KERNEL_COEFF_4 30.557749073644f
#define KERNEL_COEFF_5 5.092958178941f
#define KERNEL_COEFF_6 (-15.278874536822f)
#define NORM_COEFF 4.188790204786f
float r = sqrtf( r2 );
float ui, uj, wi, wj;
float ui_dh, uj_dh, wi_dh, wj_dh;
if ( r2 < hi*hi && pi != NULL ) {
ui = r / hi;
ui_dh = -r / hi / hi;
if ( ui < 0.5f ) {
wi = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (ui - 1.0f) * ui * ui;
wi_dh = KERNEL_COEFF_2 * ui_dh * ui * ui
+ 2 * KERNEL_COEFF_2 * (ui - 1.0f) * ui_dh * ui;
}
else {
wi = KERNEL_COEFF_5 * (1.0f - ui) * (1.0f - ui) * (1.0f - ui);
wi_dh = -3 * KERNEL_COEFF_5 * ui_dh * (1.0f - ui) * (1.0f - ui);
}
pi->count += NORM_COEFF * wi;
pi->count_dh += NORM_COEFF * wi_dh;
pi->icount += 1;
}
if ( r2 < hj*hj && pj != NULL ) {
uj = r / hj;
uj_dh = -r / hj / hj;
if ( uj < 0.5f ) {
wj = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (uj - 1.0f) * uj * uj;
wj_dh = KERNEL_COEFF_2 * uj_dh * uj * uj
+ 2 * KERNEL_COEFF_2 * (uj - 1.0f) * uj_dh * uj;
}
else {
wj = KERNEL_COEFF_5 * (1.0f - uj) * (1.0f - uj) * (1.0f - uj);
wj_dh = -3 * KERNEL_COEFF_5 * uj_dh * (1.0f - uj) * (1.0f - uj);
}
pj->count += NORM_COEFF * wj;
pj->count_dh += NORM_COEFF * wj_dh;
pj->icount += 1;
}
#ifdef HIST
if ( hi > hj )
runner_hist_hit( hi / hj );
else
runner_hist_hit( hj / hi );
#endif
}
/* A struct representing a runner's thread and its data. */ /* A struct representing a runner's thread and its data. */
struct runner_thread { struct runner_thread {
...@@ -280,7 +155,8 @@ struct runner { ...@@ -280,7 +155,8 @@ struct runner {
/* Function prototypes. */ /* Function prototypes. */
void runner_run ( struct runner *r , int sort_queues ); void runner_run ( struct runner *r , int sort_queues );
void runner_dopair ( struct runner_thread *rt , struct cell *ci , struct cell *cj ); void runner_dopair_density ( struct runner_thread *rt , struct cell *ci , struct cell *cj );
void runner_doself ( struct runner_thread *rt , struct cell *c ); void runner_doself_density ( struct runner_thread *rt , struct cell *c );
void runner_dosub_density ( struct runner_thread *rt , struct cell *ci , struct cell *cj , int flags );
void runner_dosort ( struct runner_thread *rt , struct cell *c , int flag ); void runner_dosort ( struct runner_thread *rt , struct cell *c , int flag );
void runner_init ( struct runner *r , struct space *s , int nr_threads , int nr_queues , int policy ); void runner_init ( struct runner *r , struct space *s , int nr_threads , int nr_queues , int policy );
src/runner_doiact.h 0 → 100644
+ 543
− 0
View file @ 3fd578a4
This diff is collapsed.
src/runner_iact.h 0 → 100644
+ 98
− 0
View file @ 3fd578a4
/*******************************************************************************
* This file is part of GadgetSMP.
* Coypright (c) 2012 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/>.
*
******************************************************************************/
__attribute__ ((always_inline)) INLINE void runner_iact_density ( float r2 , float hi , float hj , struct part *pi , struct part *pj ) {
#define KERNEL_COEFF_1 2.546479089470f
#define KERNEL_COEFF_2 15.278874536822f
#define KERNEL_COEFF_3 45.836623610466f
#define KERNEL_COEFF_4 30.557749073644f
#define KERNEL_COEFF_5 5.092958178941f
#define KERNEL_COEFF_6 (-15.278874536822f)
#define NORM_COEFF 4.188790204786f
float r = sqrtf( r2 );
float ui, uj, wi, wj;
float ui_dh, uj_dh, wi_dh, wj_dh;
if ( r2 < hi*hi && pi != NULL ) {
ui = r / hi;
ui_dh = -r / hi / hi;
if ( ui < 0.5f ) {
wi = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (ui - 1.0f) * ui * ui;
wi_dh = KERNEL_COEFF_2 * ui_dh * ui * ui
+ 2 * KERNEL_COEFF_2 * (ui - 1.0f) * ui_dh * ui;
}
else {
wi = KERNEL_COEFF_5 * (1.0f - ui) * (1.0f - ui) * (1.0f - ui);
wi_dh = -3 * KERNEL_COEFF_5 * ui_dh * (1.0f - ui) * (1.0f - ui);
}
pi->rho += NORM_COEFF * wi;
pi->rho_dh += NORM_COEFF * wi_dh;
pi->icount += 1;
}
if ( r2 < hj*hj && pj != NULL ) {
uj = r / hj;
uj_dh = -r / hj / hj;
if ( uj < 0.5f ) {
wj = KERNEL_COEFF_1 + KERNEL_COEFF_2 * (uj - 1.0f) * uj * uj;
wj_dh = KERNEL_COEFF_2 * uj_dh * uj * uj
+ 2 * KERNEL_COEFF_2 * (uj - 1.0f) * uj_dh * uj;
}
else {
wj = KERNEL_COEFF_5 * (1.0f - uj) * (1.0f - uj) * (1.0f - uj);
wj_dh = -3 * KERNEL_COEFF_5 * uj_dh * (1.0f - uj) * (1.0f - uj);
}
pj->rho += NORM_COEFF * wj;
pj->rho_dh += NORM_COEFF * wj_dh;
pj->icount += 1;
}
#ifdef HIST
if ( hi > hj )
runner_hist_hit( hi / hj );
else
runner_hist_hit( hj / hi );
#endif
}
__attribute__ ((always_inline)) INLINE void runner_iact_force ( float r2 , float hi , float hj , struct part *pi , struct part *pj ) {
#define KERNEL_COEFF_1 2.546479089470f
#define KERNEL_COEFF_2 15.278874536822f
#define KERNEL_COEFF_3 45.836623610466f
#define KERNEL_COEFF_4 30.557749073644f
#define KERNEL_COEFF_5 5.092958178941f
#define KERNEL_COEFF_6 (-15.278874536822f)
#define NORM_COEFF 4.188790204786f
}
src/space.c
+ 173
− 128
View file @ 3fd578a4
This diff is collapsed.
src/space.h
+ 3
− 2
View file @ 3fd578a4
...@@ -52,7 +52,7 @@ struct space { ...@@ -52,7 +52,7 @@ struct space {
double h[3], ih[3]; double h[3], ih[3];
/* The minimum and maximum cutoff radii. */ /* The minimum and maximum cutoff radii. */
double r_min, r_max; double h_min, h_max;
/* Number of cells. */ /* Number of cells. */
int nr_cells, tot_cells; int nr_cells, tot_cells;
...@@ -90,9 +90,10 @@ struct space { ...@@ -90,9 +90,10 @@ struct space {
/* function prototypes. */ /* function prototypes. */
struct cell *space_getcell ( struct space *s ); struct cell *space_getcell ( struct space *s );
struct task *space_gettask ( struct space *s ); struct task *space_gettask ( struct space *s );
struct task *space_addtask ( struct space *s , int type , int subtype , int flags , int wait , struct cell *ci , struct cell *cj , struct task *unlock_tasks[] , int nr_unlock_tasks , struct cell *unlock_cells[] , int nr_unlock_cells );
void space_init ( struct space *s , double dim[3] , struct part *parts , int N , int periodic , double h_max ); void space_init ( struct space *s , double dim[3] , struct part *parts , int N , int periodic , double h_max );
void space_maketasks ( struct space *s , int do_sort ); void space_maketasks ( struct space *s , int do_sort );
void space_map_cells ( struct space *s , void (*fun)( struct cell *c , void *data ) , void *data ); void space_map_cells ( struct space *s , int full , void (*fun)( struct cell *c , void *data ) , void *data );
void space_map_parts ( struct space *s , void (*fun)( struct part *p , struct cell *c , void *data ) , void *data ); void space_map_parts ( struct space *s , void (*fun)( struct part *p , struct cell *c , void *data ) , void *data );
void space_recycle ( struct space *s , struct cell *c ); void space_recycle ( struct space *s , struct cell *c );
......
src/task.c
+ 1
− 2
View file @ 3fd578a4
...@@ -38,8 +38,7 @@ ...@@ -38,8 +38,7 @@
/* Task type names. */ /* Task type names. */
const char *taskID_names[tid_count] = { "none" , "sort" , "self" , "pair" , "sub" }; const char *taskID_names[task_type_count] = { "none" , "sort" , "self" , "pair" , "sub" , "ghost" };
/* Error macro. */ /* Error macro. */
#define error(s) { printf( "%s:%s:%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); } #define error(s) { printf( "%s:%s:%i: %s\n" , __FILE__ , __FUNCTION__ , __LINE__ , s ); abort(); }
......
src/task.h
+ 20
− 10
View file @ 3fd578a4
...@@ -20,26 +20,36 @@ ...@@ -20,26 +20,36 @@
/* Some constants. */ /* Some constants. */
#define task_maxwait 3 #define task_maxwait 3
#define task_maxunlock 39 #define task_maxunlock 40
/* The different task IDs. */ /* The different task types. */
enum taskIDs { enum task_types {
tid_none = 0, task_type_none = 0,
tid_sort, task_type_sort,
tid_self, task_type_self,
tid_pair, task_type_pair,
tid_sub, task_type_sub,
tid_count task_type_ghost,
task_type_count
}; };
extern const char *taskID_names[]; extern const char *taskID_names[];
/* The different task sub-types. */
enum task_subtypes {
task_subtype_none = 0,
task_subtype_density,
task_subtype_force,
task_subtype_count
};
extern const char *taskID_names[];
/* Data of a task. */ /* Data of a task. */
struct task { struct task {
int type, flags, wait, rank, done; int type, subtype, flags, wait, rank, done;
int nr_unlock_tasks; int nr_unlock_tasks;
struct task *unlock_tasks[ task_maxunlock ]; struct task *unlock_tasks[ task_maxunlock ];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment