Commit 5761afd4 authored by Varan's avatar Varan
Browse files

made RANDOMIZE_DOMAINCENTER default, implemented random shifts for periodic boxes

parent 06e97192
...@@ -67,7 +67,6 @@ ...@@ -67,7 +67,6 @@
#DIRECT_SUMMATION_THRESHOLD=1000 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation #DIRECT_SUMMATION_THRESHOLD=1000 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation
#EXACT_GRAVITY_FOR_PARTICLE_TYPE=4 #N-squared fashion gravity for a small number of particles of the given type #EXACT_GRAVITY_FOR_PARTICLE_TYPE=4 #N-squared fashion gravity for a small number of particles of the given type
#EVALPOTENTIAL # computes gravitational potential #EVALPOTENTIAL # computes gravitational potential
#RANDOMIZE_DOMAINCENTER # random displacement to position of domain center; avoids correlated force-errors, important mainly for isolated systems (which otherwise might start to drift in some direction).
#--------------------------------------- TreePM Options; default: no Particle-Mesh #--------------------------------------- TreePM Options; default: no Particle-Mesh
#PMGRID=512 # Enables particle mesh; number of cells used for grid in each dimension #PMGRID=512 # Enables particle mesh; number of cells used for grid in each dimension
......
...@@ -62,6 +62,7 @@ DISABLE_SPATIAL_EXTRAPOLATION ...@@ -62,6 +62,7 @@ DISABLE_SPATIAL_EXTRAPOLATION
DISABLE_SPATIAL_RECONSTRUCTION DISABLE_SPATIAL_RECONSTRUCTION
DISABLE_TIME_EXTRAPOLATION DISABLE_TIME_EXTRAPOLATION
DISABLE_VELOCITY_CSND_SLOPE_LIMITING DISABLE_VELOCITY_CSND_SLOPE_LIMITING
DO_NOT_RANDOMIZE_DOMAINCENTER
NOCALLSOFSYSTEM NOCALLSOFSYSTEM
NOTYPEPREFIX_FFTW NOTYPEPREFIX_FFTW
NO_ID_UNIQUE_CHECK NO_ID_UNIQUE_CHECK
......
...@@ -30,7 +30,6 @@ CELL_CENTER_GRAVITY # uses geometric centers to calculate g ...@@ -30,7 +30,6 @@ CELL_CENTER_GRAVITY # uses geometric centers to calculate g
ALLOW_DIRECT_SUMMATION # Performed direct summation instead of tree-based gravity if number of active particles < DIRECT_SUMMATION_THRESHOLD (= 3000 unless specified differently here) ALLOW_DIRECT_SUMMATION # Performed direct summation instead of tree-based gravity if number of active particles < DIRECT_SUMMATION_THRESHOLD (= 3000 unless specified differently here)
DIRECT_SUMMATION_THRESHOLD=500 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation DIRECT_SUMMATION_THRESHOLD=500 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation
GRAVITY_NOT_PERIODIC # gravity is not treated periodically GRAVITY_NOT_PERIODIC # gravity is not treated periodically
RANDOMIZE_DOMAINCENTER # random displacement to position of domain center; avoids correlated force-errors, important mainly for isolated systems (which otherwise might start to drift in some direction).
#--------------------------------------- Gravity softening #--------------------------------------- Gravity softening
......
...@@ -13,7 +13,6 @@ CELL_CENTER_GRAVITY # uses geometric centers to calculate g ...@@ -13,7 +13,6 @@ CELL_CENTER_GRAVITY # uses geometric centers to calculate g
ALLOW_DIRECT_SUMMATION # Performed direct summation instead of tree-based gravity if number of active particles < DIRECT_SUMMATION_THRESHOLD (= 3000 unless specified differently here) ALLOW_DIRECT_SUMMATION # Performed direct summation instead of tree-based gravity if number of active particles < DIRECT_SUMMATION_THRESHOLD (= 3000 unless specified differently here)
DIRECT_SUMMATION_THRESHOLD=500 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation DIRECT_SUMMATION_THRESHOLD=500 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation
GRAVITY_NOT_PERIODIC # gravity is not treated periodically GRAVITY_NOT_PERIODIC # gravity is not treated periodically
RANDOMIZE_DOMAINCENTER # random displacement to position of domain center; avoids correlated force-errors, important mainly for isolated systems (which otherwise might start to drift in some direction).
#--------------------------------------- Gravity softening #--------------------------------------- Gravity softening
......
...@@ -42,6 +42,11 @@ ...@@ -42,6 +42,11 @@
#define MASK_ACTIVE_FLAG_IN_TYPE 127 #define MASK_ACTIVE_FLAG_IN_TYPE 127
#define SET_ACTIVE_FLAG_IN_TYPE 128 #define SET_ACTIVE_FLAG_IN_TYPE 128
enum domain_displace_mode
{
DISPLACE_POSITION_FORWARD,
DISPLACE_POSITION_BACKWARD
};
extern struct local_topnode_data extern struct local_topnode_data
{ {
...@@ -163,6 +168,7 @@ int domain_compare_local_trans_data_ID(const void *a, const void *b); ...@@ -163,6 +168,7 @@ int domain_compare_local_trans_data_ID(const void *a, const void *b);
int domain_compare_recv_trans_data_ID(const void *a, const void *b); int domain_compare_recv_trans_data_ID(const void *a, const void *b);
int domain_compare_recv_trans_data_oldtask(const void *a, const void *b); int domain_compare_recv_trans_data_oldtask(const void *a, const void *b);
void mysort_domain(void *b, size_t n, size_t s); void mysort_domain(void *b, size_t n, size_t s);
void domain_displacePosition( MyDouble *pos, enum domain_displace_mode mode );
#endif /* #ifndef DOMAIN_H */ #endif /* #ifndef DOMAIN_H */
...@@ -147,6 +147,14 @@ void domain_exchange_and_update_DC(void) ...@@ -147,6 +147,14 @@ void domain_exchange_and_update_DC(void)
{ {
double t0 = second(); double t0 = second();
#if !defined(GRAVITY_NOT_PERIODIC) && !defined(DO_NOT_RANDOMIZE_DOMAINCENTER) && defined(SELFGRAVITY)
/* remove all image flags, after our box movement stunt they are all incorrect anyway */
for(int i = 0; i < MaxNvc; i++)
{
DC[i].image_flags = 1;
}
#endif
/* first, we need to complete the translation table */ /* first, we need to complete the translation table */
for(int j = 0; j < NTask; j++) for(int j = 0; j < NTask; j++)
Send_count[j] = 0; Send_count[j] = 0;
...@@ -338,6 +346,9 @@ void domain_exchange_and_update_DC(void) ...@@ -338,6 +346,9 @@ void domain_exchange_and_update_DC(void)
recv_transscribe_data[i].new_task = trans_table[old_index].new_task; recv_transscribe_data[i].new_task = trans_table[old_index].new_task;
recv_transscribe_data[i].new_index = trans_table[old_index].new_index; recv_transscribe_data[i].new_index = trans_table[old_index].new_index;
#if !defined(GRAVITY_NOT_PERIODIC) && !defined(DO_NOT_RANDOMIZE_DOMAINCENTER) && defined(SELFGRAVITY)
// Nothing to do here
#else
if(recv_transscribe_data[i].new_task >= 0) if(recv_transscribe_data[i].new_task >= 0)
{ {
if(trans_table[old_index].wrapped) if(trans_table[old_index].wrapped)
...@@ -417,6 +428,7 @@ void domain_exchange_and_update_DC(void) ...@@ -417,6 +428,7 @@ void domain_exchange_and_update_DC(void)
recv_transscribe_data[i].image_flags = (1 << (zbits * 9 + ybits * 3 + xbits)); recv_transscribe_data[i].image_flags = (1 << (zbits * 9 + ybits * 3 + xbits));
} }
} }
#endif
} }
/* now return the data */ /* now return the data */
......
...@@ -45,6 +45,53 @@ ...@@ -45,6 +45,53 @@
#include "domain.h" #include "domain.h"
#include "../mesh/voronoi/voronoi.h" #include "../mesh/voronoi/voronoi.h"
/*! \brief Move the coordinate in pos by the global displacement vector
*
* \param[in] pos coordinate vector (3 entries).
* \param[in] mode displacement mode, either DISPLACE_POSITION_FORWARD or DISPLACE_POSITION_BACKWARD
*
* \return void
*/
void domain_displacePosition( MyDouble *pos, enum domain_displace_mode mode )
{
if(mode == DISPLACE_POSITION_FORWARD)
{
double xtmp, ytmp, ztmp;
pos[0] = WRAP_X( pos[0] + All.GlobalDisplacementVector[0] );
pos[1] = WRAP_Y( pos[1] + All.GlobalDisplacementVector[1] );
pos[2] = WRAP_Z( pos[2] + All.GlobalDisplacementVector[2] );
}
else if(mode == DISPLACE_POSITION_BACKWARD)
{
double xtmp, ytmp, ztmp;
pos[0] = WRAP_X( pos[0] - All.GlobalDisplacementVector[0] );
pos[1] = WRAP_Y( pos[1] - All.GlobalDisplacementVector[1] );
pos[2] = WRAP_Z( pos[2] - All.GlobalDisplacementVector[2] );
}
else
terminate( "Unkown mode %d.", mode );
}
/*! \brief Move the coordinate for all positions by the global displacement vector
*
* \param[in] mode displacement mode, either DISPLACE_POSITION_FORWARD or DISPLACE_POSITION_BACKWARD
*
* \return void
*/
static void domain_displacePositions( enum domain_displace_mode mode )
{
for(int i = 0; i < NumPart; i++)
{
if(P[i].ID == 0 && P[i].Mass == 0) /* derefined */
continue;
domain_displacePosition( P[i].Pos, mode );
if(i < NumGas)
domain_displacePosition( SphP[i].Center, mode );
}
}
/*! \brief Finds the extent of the global domain grid. /*! \brief Finds the extent of the global domain grid.
* *
...@@ -116,7 +163,7 @@ void domain_findExtent(void) ...@@ -116,7 +163,7 @@ void domain_findExtent(void)
#endif /* #if defined(GRAVITY_NOT_PERIODIC) && !defined(ADDBACKGROUNDGRID) #else */ #endif /* #if defined(GRAVITY_NOT_PERIODIC) && !defined(ADDBACKGROUNDGRID) #else */
#if !defined(RANDOMIZE_DOMAINCENTER) #if defined(DO_NOT_RANDOMIZE_DOMAINCENTER) || !defined(GRAVITY_NOT_PERIODIC) || defined(ONEDIMS) || defined(TWODIMS)
for(j = 0; j < 3; j++) for(j = 0; j < 3; j++)
{ {
DomainCenter[j] = 0.5 * (xmin_glob[j] + xmax_glob[j]); DomainCenter[j] = 0.5 * (xmin_glob[j] + xmax_glob[j]);
...@@ -176,6 +223,19 @@ void do_box_wrapping(void) ...@@ -176,6 +223,19 @@ void do_box_wrapping(void)
boxsize[2] *= LONG_Z; boxsize[2] *= LONG_Z;
#endif /* #ifdef LONG_Z */ #endif /* #ifdef LONG_Z */
#if !defined(GRAVITY_NOT_PERIODIC) && !defined(DO_NOT_RANDOMIZE_DOMAINCENTER) && defined(SELFGRAVITY) && (NUMDIMS > 2)
domain_displacePositions( DISPLACE_POSITION_BACKWARD );
if(ThisTask == 0)
for(j = 0; j < 3; j++)
All.GlobalDisplacementVector[j] = (get_random_number() - 0.5) * boxsize[j];
mpi_printf("DOMAIN: New global displacement vector: %g, %g, %g\n", All.GlobalDisplacementVector[0], All.GlobalDisplacementVector[1], All.GlobalDisplacementVector[2]);
MPI_Bcast(All.GlobalDisplacementVector, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
domain_displacePositions( DISPLACE_POSITION_FORWARD );
#endif
int i; int i;
for(i = 0; i < NumPart; i++) for(i = 0; i < NumPart; i++)
{ {
......
...@@ -83,6 +83,53 @@ void write_compile_time_options_in_hdf5(hid_t handle); ...@@ -83,6 +83,53 @@ void write_compile_time_options_in_hdf5(hid_t handle);
#ifdef FOF #ifdef FOF
/*! \brief Make sure a position lies in the box in case of periodic boundaries.
*
* \param[in] pos Single coordinate in one dimension to be wrapped
* \param[in] dim Index of coordinate [0/1/2]
*
* \return double: wrapped coordinate
*/
MyOutputFloat static wrap_position( MyOutputFloat pos, int dim )
{
#if defined(REFLECTIVE_X)
if(dim == 0)
return pos;
#endif
#if defined(REFLECTIVE_Y)
if(dim == 1)
return pos;
#endif
#if defined(REFLECTIVE_Z)
if(dim == 2)
return pos;
#endif
double boxsize = All.BoxSize;
#ifdef LONG_X
if(dim == 0)
boxsize *= LONG_X;
#endif
#ifdef LONG_Y
if(dim == 1)
boxsize *= LONG_Y;
#endif
#ifdef LONG_Z
if(dim == 2)
boxsize *= LONG_Z;
#endif
while(pos < 0)
pos += boxsize;
while(pos >= boxsize)
pos -= boxsize;
return pos;
}
/*! \brief Main routine for group output. /*! \brief Main routine for group output.
* *
...@@ -576,14 +623,14 @@ void fof_subfind_fill_write_buffer(enum fof_subfind_iofields blocknr, int *start ...@@ -576,14 +623,14 @@ void fof_subfind_fill_write_buffer(enum fof_subfind_iofields blocknr, int *start
case IO_FOF_POS: case IO_FOF_POS:
for(k = 0; k < 3; k++) for(k = 0; k < 3; k++)
#ifdef SUBFIND #ifdef SUBFIND
*fp++ = Group[pindex].Pos[k]; *fp++ = wrap_position( Group[pindex].Pos[k] - All.GlobalDisplacementVector[k], k );
#else /* #ifdef SUBFIND */ #else /* #ifdef SUBFIND */
*fp++ = Group[pindex].CM[k]; *fp++ = wrap_position( Group[pindex].CM[k] - All.GlobalDisplacementVector[k], k );
#endif /* #ifdef SUBFIND #else */ #endif /* #ifdef SUBFIND #else */
break; break;
case IO_FOF_CM: case IO_FOF_CM:
for(k = 0; k < 3; k++) for(k = 0; k < 3; k++)
*fp++ = Group[pindex].CM[k]; *fp++ = wrap_position( Group[pindex].CM[k] - All.GlobalDisplacementVector[k], k );
break; break;
case IO_FOF_VEL: case IO_FOF_VEL:
for(k = 0; k < 3; k++) for(k = 0; k < 3; k++)
...@@ -1094,7 +1141,7 @@ void fof_subfind_fill_write_buffer(enum fof_subfind_iofields blocknr, int *start ...@@ -1094,7 +1141,7 @@ void fof_subfind_fill_write_buffer(enum fof_subfind_iofields blocknr, int *start
case IO_SUB_POS: case IO_SUB_POS:
#ifdef SUBFIND #ifdef SUBFIND
for(k = 0; k < 3; k++) for(k = 0; k < 3; k++)
*fp++ = SubGroup[pindex].Pos[k]; *fp++ = wrap_position( SubGroup[pindex].Pos[k] - All.GlobalDisplacementVector[k], k );
#endif /* #ifdef SUBFIND */ #endif /* #ifdef SUBFIND */
break; break;
case IO_SUB_VEL: case IO_SUB_VEL:
...@@ -1118,7 +1165,7 @@ void fof_subfind_fill_write_buffer(enum fof_subfind_iofields blocknr, int *start ...@@ -1118,7 +1165,7 @@ void fof_subfind_fill_write_buffer(enum fof_subfind_iofields blocknr, int *start
case IO_SUB_CM: case IO_SUB_CM:
#ifdef SUBFIND #ifdef SUBFIND
for(k = 0; k < 3; k++) for(k = 0; k < 3; k++)
*fp++ = SubGroup[pindex].CM[k]; *fp++ = wrap_position( SubGroup[pindex].CM[k] - All.GlobalDisplacementVector[k], k );
#endif /* #ifdef SUBFIND */ #endif /* #ifdef SUBFIND */
break; break;
case IO_SUB_SPIN: case IO_SUB_SPIN:
......
...@@ -92,6 +92,9 @@ int init(void) ...@@ -92,6 +92,9 @@ int init(void)
} }
set_cosmo_factors_for_current_time(); set_cosmo_factors_for_current_time();
for(j = 0; j < 3; j++)
All.GlobalDisplacementVector[j] = 0;
All.NumCurrentTiStep = 0; /* setup some counters */ All.NumCurrentTiStep = 0; /* setup some counters */
All.SnapshotFileCount = 0; All.SnapshotFileCount = 0;
......
...@@ -170,7 +170,7 @@ void io_func_pos(int particle, int components, void *buffer, int mode) ...@@ -170,7 +170,7 @@ void io_func_pos(int particle, int components, void *buffer, int mode)
for(k = 0; k < 3; k++) for(k = 0; k < 3; k++)
{ {
pp[k] = P[particle].Pos[k]; pp[k] = P[particle].Pos[k] - All.GlobalDisplacementVector[k];
#if defined(GRAVITY_NOT_PERIODIC) #if defined(GRAVITY_NOT_PERIODIC)
if(P[particle].Type != 0) if(P[particle].Type != 0)
...@@ -202,7 +202,7 @@ void io_func_pos(int particle, int components, void *buffer, int mode) ...@@ -202,7 +202,7 @@ void io_func_pos(int particle, int components, void *buffer, int mode)
for(k = 0; k < 3; k++) for(k = 0; k < 3; k++)
{ {
pp[k] = P[particle].Pos[k]; pp[k] = P[particle].Pos[k] - All.GlobalDisplacementVector[k];
#if defined(GRAVITY_NOT_PERIODIC) #if defined(GRAVITY_NOT_PERIODIC)
if(P[particle].Type != 0) if(P[particle].Type != 0)
...@@ -238,7 +238,7 @@ void io_func_pos(int particle, int components, void *buffer, int mode) ...@@ -238,7 +238,7 @@ void io_func_pos(int particle, int components, void *buffer, int mode)
for(k = 0; k < components; k++) for(k = 0; k < components; k++)
{ {
P[particle].Pos[k] = in_buffer[k]; P[particle].Pos[k] = in_buffer[k] + All.GlobalDisplacementVector[k];
} }
} }
} }
......
...@@ -1287,6 +1287,7 @@ MyIDType MaxID; ...@@ -1287,6 +1287,7 @@ MyIDType MaxID;
double CoreRadius; double CoreRadius;
#endif /* #ifdef ONEDIMS_SPHERICAL */ #endif /* #ifdef ONEDIMS_SPHERICAL */
double GlobalDisplacementVector[3];
} }
All; All;
......
Markdown is supported
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