6 Commits

Author SHA1 Message Date
be0e67e954 using momentum changes now works. 2019-07-16 21:22:11 +02:00
26acc9802d fixed a typo 2019-07-16 21:19:04 +02:00
1e9bc3209c merged origin/use_momentum_changes 2019-07-16 21:02:53 +02:00
fd6e4d6bc4 some work 2019-07-16 21:02:17 +02:00
9580818204 merged 2019-07-15 22:31:53 +02:00
1a9ddddc6d some work; does not work 2019-07-15 16:56:16 +02:00
2 changed files with 17 additions and 21 deletions

View File

@@ -113,17 +113,13 @@ interaction_ufunc_float2D
char * x_old = args[0] char * x_old = args[0]
, * y_old = args[1] , * y_old = args[1]
, * p_x_old = args[2] , * p_x_new = args[2]
, * p_y_old = args[3] , * p_y_new = args[3];
, * p_x_new = args[4]
, * p_y_new = args[5];
npy_intp x_old_steps = steps[0] npy_intp x_old_steps = steps[0]
, y_old_steps = steps[1] , y_old_steps = steps[1]
, p_x_old_steps = steps[2] , p_x_new_steps = steps[2]
, p_y_old_steps = steps[3] , p_y_new_steps = steps[3];
, p_x_new_steps = steps[4]
, p_y_new_steps = steps[5];
float * coefficients = (float *) data; float * coefficients = (float *) data;
float dt = coefficients[19]; float dt = coefficients[19];
@@ -133,7 +129,7 @@ interaction_ufunc_float2D
// Stuff we will need: // Stuff we will need:
float r; // Distance between interacting particles. float r; // Distance between interacting particles.
float delta_x_e; // x-component of unit direction vector. float delta_x_e; // x-component of unit direction vector.
float delty_y_e; // y-component of unit direction vector. float delta_y_e; // y-component of unit direction vector.
float delta_p_x; float delta_p_x;
float delta_p_y; float delta_p_y;
@@ -151,10 +147,8 @@ interaction_ufunc_float2D
this_x_i = *(float *)(x_old + i*x_old_steps); this_x_i = *(float *)(x_old + i*x_old_steps);
this_y_i = *(float *)(y_old + i*y_old_steps); this_y_i = *(float *)(y_old + i*y_old_steps);
// copy current momenta *(float *)(p_x_new + i*p_x_new_steps) = 0;
*(float *)(p_x_new + i*p_x_new_steps) = *(float *)(p_x_old + i*p_x_old_steps); *(float *)(p_y_new + i*p_y_new_steps) = 0;
*(float *)(p_y_new + i*p_y_new_steps) = *(float *)(p_y_old + i*p_y_old_steps);
// compute and add the momentum offset // compute and add the momentum offset
for(j = 0; j < i; j++) for(j = 0; j < i; j++)
{ {
@@ -168,7 +162,7 @@ interaction_ufunc_float2D
if(r > 0) if(r > 0)
{ {
delta_x_e = (this_x_i - this_x_j) / r; delta_x_e = (this_x_i - this_x_j) / r;
delty_y_e = (this_y_i - this_y_j) / r; delta_y_e = (this_y_i - this_y_j) / r;
} }
else else
{ {
@@ -178,13 +172,13 @@ interaction_ufunc_float2D
printf("Warning: corrected r = 0 with random angle pi*%f\n", random_angle / M_PI); printf("Warning: corrected r = 0 with random angle pi*%f\n", random_angle / M_PI);
#endif #endif
delta_x_e = cosf(random_angle); delta_x_e = cosf(random_angle);
delty_y_e = sinf(random_angle); delta_y_e = sinf(random_angle);
} }
// Update the momenta. // Update the momenta.
delta_p_x = delta_x_e * interaction_force_function(r, coefficients); delta_p_x = delta_x_e * interaction_force_function(r, coefficients);
delta_p_y = delty_y_e * interaction_force_function(r, coefficients); delta_p_y = delta_y_e * interaction_force_function(r, coefficients);
#ifdef WARN_NAN_OCCUR #ifdef WARN_NAN_OCCUR
if(isnan(delta_p_x)) if(isnan(delta_p_x))
{ {
@@ -205,7 +199,7 @@ interaction_ufunc_float2D
//NPY_END_THREADS; //NPY_END_THREADS;
} }
static char interaction_types[] = static char interaction_types[] =
{ NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT}; { NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT};
static char force_types[] = static char force_types[] =
{ NPY_FLOAT, NPY_FLOAT}; { NPY_FLOAT, NPY_FLOAT};
static char potential_types[] = static char potential_types[] =
@@ -298,9 +292,9 @@ interaction_UFuncWrapper_init
interaction_funcs interaction_funcs
, self->data , self->data
, interaction_types , interaction_types
, 1 , 1 // ntypes
, 4 , 2 // nin
, 2 , 2 // nout
, PyUFunc_None , PyUFunc_None
, "interaction2D" , "interaction2D"
, "Update the momenta according to the given coefficients and positions" , "Update the momenta according to the given coefficients and positions"

View File

@@ -37,7 +37,9 @@ class BrownIterator(object):
if(self._i == 1): if(self._i == 1):
return self.x, self.y return self.x, self.y
self.px, self.py = self._interaction(self.x, self.y, self.px, self.py) delta_px, delta_py = self._interaction(self.x, self.y)
self.px = delta_px + self.px
self.py = delta_py + self.py
# XXX: We need the (-1)**i to make the problem # XXX: We need the (-1)**i to make the problem
# symmetric. # symmetric.
self.px[np.isnan(self.px)] = self.speed_of_light * (-1)**self._i self.px[np.isnan(self.px)] = self.speed_of_light * (-1)**self._i