10 Commits

Author SHA1 Message Date
a9450179be Merge branch 'master' into use_momentum_changes 2019-07-16 21:22:22 +02:00
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
609c9c3d1a some nice potential 2019-07-16 21:16:08 +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
48b8ce8c88 some more useful comments 2019-07-16 20:45:17 +02:00
da2fd89946 added some useful comments 2019-07-16 20:44:15 +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
5 changed files with 41 additions and 23 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

@@ -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,9 +4,12 @@ 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)
# Plot the force and potential.
r = np.arange(0, 100, 0.02, dtype=np.float16) r = np.arange(0, 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")

View File

@@ -11,23 +11,34 @@ 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 = 600
# Idk, seems to not do anyting.
frames = 100 frames = 100
# Only spawn in 1/x of the borders.
spawn_restriction = 3 spawn_restriction = 3
# 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.1 dt = 0.1
c[-1] = dt c[-1] = dt
# Initial positions.
x_coords = np.random.uniform(borders_x[0] / spawn_restriction, borders_x[1] / spawn_restriction, n_particles).astype(np.float16) x_coords = np.random.uniform(borders_x[0] / spawn_restriction, borders_x[1] / spawn_restriction, n_particles).astype(np.float16)
y_coords = np.random.uniform(borders_y[0] / spawn_restriction, borders_y[1] / spawn_restriction, n_particles).astype(np.float16) y_coords = np.random.uniform(borders_y[0] / spawn_restriction, borders_y[1] / spawn_restriction, n_particles).astype(np.float16)
# Initial momenta are 0.
x_momenta = np.zeros(n_particles, dtype=np.float16) x_momenta = np.zeros(n_particles, dtype=np.float16)
y_momenta = np.zeros(n_particles, dtype=np.float16) y_momenta = np.zeros(n_particles, dtype=np.float16)
# 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,22 +46,28 @@ 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-")
# 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):
# 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())

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