From d2336cded2967877b74c2a9c6ee815f8d470b65d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Mon, 15 Jul 2019 21:45:32 +0200 Subject: [PATCH 1/2] initial changes towards using the potential --- c/interaction/interaction.c | 65 +++++++++++++++++++++++++++++++++++++ coefficients.py | 8 ++--- force.py | 5 ++- particles.py | 4 +-- 4 files changed, 75 insertions(+), 7 deletions(-) diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index 7da2163..41b0728 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -10,6 +10,28 @@ static float interaction_force_function ( float r , float * coefficients) +{ + float result = 0; + int i; + result = coefficients[0]; + for(i = 1; i < 7; i++) + { + result += coefficients[i] * (powf(r, i - 1) + coefficients[8]*powf(r, i)); + } + result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9])); + result += coefficients[10] * coefficients[11] * expf(coefficients[11] * (r - coefficients[12])); + result += coefficients[13] + * 2 * (r - coefficients[15]) + * coefficients[14] * expf(coefficients[14] * raise2(r - coefficients[15])); + result += coefficients[16] + * 2 * (r - coefficients[18]) + * coefficients[17] * expf(coefficients[17] * raise2(r - coefficients[18])); + return result; +} +static float +interaction_potential_function + ( float r + , float * coefficients) { float result = 0; int i; @@ -49,6 +71,30 @@ interaction_ufunc_force } } +static void +interaction_ufunc_potential + ( char ** args + , npy_intp * dimensions + , npy_intp * steps + , void * data) +{ + char * in = args[0] + , * out = args[1]; + npy_intp n = dimensions[0]; + npy_intp in_step = steps[0] + , out_step = steps[1]; + + npy_intp i; + + float * coefficients = (float *) data; + + for(i = 0; i < n; i++) + { + *(float *)out = interaction_potential_function(*(float *)in, coefficients); + out += out_step; + in += in_step; + } +} static void interaction_ufunc_float2D @@ -159,6 +205,10 @@ static char interaction_types[] = { NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT, NPY_FLOAT}; static char force_types[] = { NPY_FLOAT, NPY_FLOAT}; +static char potential_types[] = + { NPY_FLOAT, NPY_FLOAT}; +static PyUFuncGenericFunction potential_funcs[1] = + { interaction_ufunc_potential}; static PyUFuncGenericFunction force_funcs[1] = { interaction_ufunc_force}; static PyUFuncGenericFunction interaction_funcs[1] = @@ -254,6 +304,21 @@ interaction_UFuncWrapper_init , 0); break; } + case 2: + { + self->ufunc = PyUFunc_FromFuncAndData( + potential_funcs // func + , self->data // data + , potential_types //types + , 1 // ntypes + , 1 // nin + , 1 // nout + , PyUFunc_None // identity + , "potential_function" // name + , "computes the scalar potential between two particles with given coefficients" // doc + , 0); // unused + break; + } default: { PyErr_SetString(PyExc_ValueError, "unknown ufunc type, must be 0 or 1"); diff --git a/coefficients.py b/coefficients.py index e7bac49..0056c67 100644 --- a/coefficients.py +++ b/coefficients.py @@ -8,15 +8,15 @@ c = np.array( , 0 # x**5 , 0 # x**6) , 1 # *c*exp( - , -1 # c + , -0.7 # c , 0 # (r - c)) - , -1 # + c*exp( + , 0 # + c*exp( , -.1 # c , 2 # (r - c)) - , -4 # + c*exp( + , -0 # + c*exp( , -0.05 # c , 40 # (r - c)**2) - , 10 # + c*exp( + , 0 # + c*exp( , -0.05 # c , 20 # (r - c)**2) , 1] # dt diff --git a/force.py b/force.py index 97eedc4..e6742c1 100644 --- a/force.py +++ b/force.py @@ -5,7 +5,10 @@ import matplotlib.pyplot as plt from coefficients import c force_function = UFuncWrapper(0, c) +potential_function = UFuncWrapper(2, c) r = np.arange(0, 100, 0.02, dtype=np.float16) -plt.plot(r, force_function(r)) +f, = plt.plot(r, force_function(r), label="force") +p, = plt.plot(r, potential_function(r), label="potential") +plt.legend(handles=[f, p]) plt.show() diff --git a/particles.py b/particles.py index 7fbfe01..7718e0a 100644 --- a/particles.py +++ b/particles.py @@ -15,8 +15,8 @@ borders_x = [-100, 100] borders_y = [-100, 100] n_particles = 600 frames = 100 -spawn_restriction = 1.1 -dt = 0.01 +spawn_restriction = 3 +dt = 0.1 c[-1] = dt x_coords = np.random.uniform(borders_x[0] / spawn_restriction, borders_x[1] / spawn_restriction, n_particles).astype(np.float16) From 5a6e039e288ba84293332164e7d5f513fb6f02e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Kn=C3=BCttel?= Date: Mon, 15 Jul 2019 22:27:12 +0200 Subject: [PATCH 2/2] potentials work --- c/interaction/interaction.c | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/c/interaction/interaction.c b/c/interaction/interaction.c index 41b0728..7d11dce 100644 --- a/c/interaction/interaction.c +++ b/c/interaction/interaction.c @@ -13,19 +13,22 @@ interaction_force_function { float result = 0; int i; - result = coefficients[0]; + result = coefficients[0] * coefficients[8]; for(i = 1; i < 7; i++) { - result += coefficients[i] * (powf(r, i - 1) + coefficients[8]*powf(r, i)); + result += coefficients[i] * (powf(r, i - 1) / i + coefficients[8]*powf(r, i)); } result *= coefficients[7] * expf(coefficients[8] * (r - coefficients[9])); - result += coefficients[10] * coefficients[11] * expf(coefficients[11] * (r - coefficients[12])); + result += coefficients[10] * coefficients[11] + * expf(coefficients[11] * (r - coefficients[12])); result += coefficients[13] * 2 * (r - coefficients[15]) - * coefficients[14] * expf(coefficients[14] * raise2(r - coefficients[15])); + * coefficients[14] + * expf(coefficients[14] * raise2(r - coefficients[15])); result += coefficients[16] * 2 * (r - coefficients[18]) - * coefficients[17] * expf(coefficients[17] * raise2(r - coefficients[18])); + * coefficients[17] + * expf(coefficients[17] * raise2(r - coefficients[18])); return result; } static float @@ -193,10 +196,10 @@ interaction_ufunc_float2D } #endif //printf("%d, %d: %f, %f\n", i, j, delta_p_x, delta_p_y); - *(float *)(p_x_new + i*p_x_new_steps) += dt*delta_p_x; - *(float *)(p_y_new + i*p_y_new_steps) += dt*delta_p_y; - *(float *)(p_x_new + j*p_x_new_steps) -= dt*delta_p_x; - *(float *)(p_y_new + j*p_y_new_steps) -= dt*delta_p_y; + *(float *)(p_x_new + i*p_x_new_steps) -= dt*delta_p_x; + *(float *)(p_y_new + i*p_y_new_steps) -= dt*delta_p_y; + *(float *)(p_x_new + j*p_x_new_steps) += dt*delta_p_x; + *(float *)(p_y_new + j*p_y_new_steps) += dt*delta_p_y; } } //NPY_END_THREADS;