Compare commits

..

18 Commits

7 changed files with 111 additions and 47 deletions

View File

@ -1,11 +1,10 @@
#include "interaction.h" #include "interaction.h"
#include "math.h" #include <math.h>
//#define WARN_NAN_OCCUR //#define WARN_NAN_OCCUR
//#define WARN_CORRECTED_R_0 //#define WARN_CORRECTED_R_0
#define raise2(x) (x)*(x) #define raise2(x) ((x)*(x))
static float static float
interaction_force_function interaction_force_function
( float r ( float r
@ -29,6 +28,12 @@ interaction_force_function
* 2 * (r - coefficients[18]) * 2 * (r - coefficients[18])
* coefficients[17] * coefficients[17]
* expf(coefficients[17] * raise2(r - coefficients[18])); * expf(coefficients[17] * raise2(r - coefficients[18]));
// 1 / r^2 migth produce NaN, avoid that if the term is disabled
// anyways.
if(coefficients[19])
{
result += -2 * coefficients[19] / powf(r + coefficients[20], 3);
}
return result; return result;
} }
static float static float
@ -46,6 +51,12 @@ interaction_potential_function
result += coefficients[10] * expf(coefficients[11] * (r - coefficients[12])); result += coefficients[10] * expf(coefficients[11] * (r - coefficients[12]));
result += coefficients[13] * expf(coefficients[14] * raise2(r - coefficients[15])); result += coefficients[13] * expf(coefficients[14] * raise2(r - coefficients[15]));
result += coefficients[16] * expf(coefficients[17] * raise2(r - coefficients[18])); result += coefficients[16] * expf(coefficients[17] * raise2(r - coefficients[18]));
// 1 / r migth produce NaN, avoid that if the term is disabled
// anyways.
if(coefficients[19])
{
result += coefficients[19] / raise2(r + coefficients[20]);
}
return result; return result;
} }
@ -99,6 +110,7 @@ interaction_ufunc_potential
} }
} }
static void static void
interaction_ufunc_float2D interaction_ufunc_float2D
( char ** args ( char ** args
@ -113,27 +125,23 @@ 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[21];
// Compute the new momenta: // Compute the new momenta:
// 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 +159,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 +174,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 +184,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 +211,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[] =
@ -220,7 +226,7 @@ static PyUFuncGenericFunction interaction_funcs[1] =
typedef struct typedef struct
{ {
PyObject_HEAD PyObject_HEAD
float coefficients[20]; float coefficients[22];
PyObject * ufunc; PyObject * ufunc;
void *data[1]; void *data[1];
} interaction_UFuncWrapper; } interaction_UFuncWrapper;
@ -249,14 +255,14 @@ interaction_UFuncWrapper_init
return -1; return -1;
} }
if(PySequence_Size(coefficients) != 20) if(PySequence_Size(coefficients) != 22)
{ {
PyErr_SetString(PyExc_ValueError, "coefficients must have length 20"); PyErr_SetString(PyExc_ValueError, "coefficients must have length 22");
return -1; return -1;
} }
// copy the coefficients. // copy the coefficients.
for(i = 0; i < 20; i++) for(i = 0; i < 22; i++)
{ {
this_coefficient = PySequence_GetItem(coefficients, i); this_coefficient = PySequence_GetItem(coefficients, i);
if(!this_coefficient) if(!this_coefficient)
@ -298,9 +304,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

@ -9,6 +9,21 @@
#include <numpy/ufuncobject.h> #include <numpy/ufuncobject.h>
#include <stddef.h> #include <stddef.h>
static void
interaction_ufunc_float2D
( char ** args
, npy_intp * dimensions
, npy_intp * steps
, void * data);
static void
interaction_ufunc_force
( char ** args
, npy_intp * dimensions
, npy_intp * steps
, void * data);
/* /*
* This is a quite generic force function mapping a * This is a quite generic force function mapping a
* distance to the magnitude of a force. The coefficients * distance to the magnitude of a force. The coefficients
@ -25,12 +40,10 @@ interaction_force_function
( float r ( float r
, float * coefficients); , float * coefficients);
static void static float
interaction_ufunc_float2D interaction_potential_function
( char ** args ( float r
, npy_intp * dimensions , float * coefficients);
, npy_intp * steps
, void * data);
static void static void
interaction_ufunc_force interaction_ufunc_force
@ -38,5 +51,10 @@ interaction_ufunc_force
, npy_intp * dimensions , npy_intp * dimensions
, npy_intp * steps , npy_intp * steps
, void * data); , void * data);
static void
interaction_ufunc_potential
( char ** args
, npy_intp * dimensions
, npy_intp * steps
, void * data);
#endif #endif

View File

@ -10,7 +10,7 @@ c = np.array(
, 1 # *c*exp( , 1 # *c*exp(
, -0.7 # c , -0.7 # c
, 0 # (r - c)) , 0 # (r - c))
, 0 # + c*exp( , +5 # + c*exp(
, -.1 # c , -.1 # c
, 2 # (r - c)) , 2 # (r - c))
, -0 # + c*exp( , -0 # + c*exp(

View File

@ -4,10 +4,13 @@ import matplotlib.pyplot as plt
from coefficients import c from coefficients import c
# This is the quite unreadable way to create
# UFuncs with given parameters. FIXME: add this to another module.
force_function = UFuncWrapper(0, c) force_function = UFuncWrapper(0, c)
potential_function = UFuncWrapper(2, c) potential_function = UFuncWrapper(2, c)
r = np.arange(0, 100, 0.02, dtype=np.float16) # Plot the force and potential.
r = np.arange(0.01, 100, 0.02, dtype=np.float16)
f, = plt.plot(r, force_function(r), label="force") f, = plt.plot(r, force_function(r), label="force")
p, = plt.plot(r, potential_function(r), label="potential") p, = plt.plot(r, potential_function(r), label="potential")
plt.legend(handles=[f, p]) plt.legend(handles=[f, p])

View File

@ -11,23 +11,37 @@ from coefficients import c
#force_function = UFuncWrapper(0, c) #force_function = UFuncWrapper(0, c)
#interaction2D = UFuncWrapper(1, c) #interaction2D = UFuncWrapper(1, c)
# Borders for both the plot and the boundary condition
# (the boundary condition might be deactivated, when creating
# the BrownIterator).
borders_x = [-100, 100] borders_x = [-100, 100]
borders_y = [-100, 100] borders_y = [-100, 100]
n_particles = 600 n_particles = 60
# Idk, seems to not do anyting.
frames = 100 frames = 100
spawn_restriction = 3 # Only spawn in 1/x of the borders.
dt = 0.1 spawn_restriction = 2
# Time resolution. Note that setting this to a too
# high value (i.e. low resolution) will lead to
# erratic behaviour, because potentials can be skipped.
dt = 0.01
c[-1] = dt c[-1] = dt
x_coords = np.random.uniform(borders_x[0] / spawn_restriction, borders_x[1] / spawn_restriction, n_particles).astype(np.float16) # Draw arrows, makes the simulation fucking slow.
y_coords = np.random.uniform(borders_y[0] / spawn_restriction, borders_y[1] / spawn_restriction, n_particles).astype(np.float16) draw_arrows = True
# Initial positions.
x_coords = np.random.uniform(borders_x[0] / spawn_restriction, borders_x[1] / spawn_restriction, n_particles).astype(np.float32)
y_coords = np.random.uniform(borders_y[0] / spawn_restriction, borders_y[1] / spawn_restriction, n_particles).astype(np.float32)
x_momenta = np.zeros(n_particles, dtype=np.float16) # Initial momenta are 0.
y_momenta = np.zeros(n_particles, dtype=np.float16) x_momenta = np.zeros(n_particles, dtype=np.float32)
y_momenta = np.zeros(n_particles, dtype=np.float32)
# Prepare the plot, remove axis & stuff.
fig = plt.figure(figsize=(7, 7)) fig = plt.figure(figsize=(7, 7))
ax = fig.add_axes([0, 0, 1, 1], frameon=False) ax = fig.add_axes([0, 0, 1, 1], frameon=False)
ax.set_xlim(*borders_x) ax.set_xlim(*borders_x)
@ -35,28 +49,47 @@ ax.set_xticks([])
ax.set_ylim(*borders_y) ax.set_ylim(*borders_y)
ax.set_yticks([]) ax.set_yticks([])
# Plot the initial values.
plot, = ax.plot(x_coords, y_coords, "b.") plot, = ax.plot(x_coords, y_coords, "b.")
center_of_mass, = ax.plot(x_coords.mean(), y_coords.mean(), "r-") center_of_mass, = ax.plot(x_coords.mean(), y_coords.mean(), "r-")
if(draw_arrows):
arrows = [ plt.Arrow(x, y, dx, dy) for x, y, dx, dy in zip(x_coords, y_coords, x_momenta, y_momenta)]
for arrow in arrows:
ax.add_patch(arrow)
# Keep track of the center of mass.
center_of_mass_history_x = deque([x_coords.mean()]) center_of_mass_history_x = deque([x_coords.mean()])
center_of_mass_history_y = deque([y_coords.mean()]) center_of_mass_history_y = deque([y_coords.mean()])
brown = BrownIterator(-1, c brown = BrownIterator(-1, c # Max iterations, simulation parameters.
, x_coords, y_coords , x_coords, y_coords
, y_momenta, y_momenta , y_momenta, y_momenta
# The boundary condition: reflect at the borders,
, borders_x, borders_y , borders_x, borders_y
# or just let propagate to infinity.
#, [], [] #, [], []
# Let the border dampen the system, border_dampening < 1 => energy is absorbed.
, border_dampening=1 , border_dampening=1
, dt=dt) , dt=dt)
u = iter(brown) u = iter(brown)
def update(i): def update(i):
global arrows
# Get the next set of positions.
data = next(u) data = next(u)
center_of_mass_history_x.append(x_coords.mean()) center_of_mass_history_x.append(x_coords.mean())
center_of_mass_history_y.append(y_coords.mean()) center_of_mass_history_y.append(y_coords.mean())
plot.set_data(*data) plot.set_data(*data)
center_of_mass.set_data(center_of_mass_history_x, center_of_mass_history_y) center_of_mass.set_data(center_of_mass_history_x, center_of_mass_history_y)
if(draw_arrows):
for arrow in arrows:
ax.patches.remove(arrow)
arrows = [plt.Arrow(x, y, dx, dy) for x, y, dx, dy in zip(data[0], data[1], u.px, u.py)]
for arrow in arrows:
ax.add_patch(arrow)
animation = ani.FuncAnimation(fig, update, range(frames), interval=1) animation = ani.FuncAnimation(fig, update, range(frames), interval=1)
plt.show() plt.show()

View File

@ -37,11 +37,15 @@ 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.
# FIXME: is this necessary?
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
self.py[np.isnan(self.py)] = self.speed_of_light * (-1)**self._i self.py[np.isnan(self.py)] = self.speed_of_light * (-1)**self._i
self._reflect_at_borders() self._reflect_at_borders()
self.x += self.dt * self.px self.x += self.dt * self.px

View File

@ -1,4 +1,4 @@
from distutils.core import setup, Extension from setuptools import setup,Extension
interaction = Extension("brown.interaction", interaction = Extension("brown.interaction",