Skip to content

Commit

Permalink
Addition of 2D vortex particle methods
Browse files Browse the repository at this point in the history
Includes refactoring to avoid name collisions with 3D vortex particle methods,
and updates to method names to better reflect new naming schemes.
  • Loading branch information
hjabird committed Jun 11, 2019
1 parent 20a9429 commit 2017561
Show file tree
Hide file tree
Showing 9 changed files with 708 additions and 123 deletions.
61 changes: 53 additions & 8 deletions include/cvortex/libcvtx.h
Original file line number Diff line number Diff line change
Expand Up @@ -48,15 +48,38 @@ typedef struct {
float strength; /* Vort per unit length */
} cvtx_F3D;

/* A Vortex particle/filament in 2D */
typedef struct {
float(*g_fn)(float rho);
float(*zeta_fn)(float rho);
void(*combined_fn)(float rho, float* g, float* zeta);
float(*eta_fn)(float rho);
bsv_V2f coord;
float vorticity;
float area;
} cvtx_P2D;

/* Vortex particle regularisation functions
Naming is following that of Winckelmans
- g(rho): normally used in induced vel
(NULL for unsupported)
- zeta(rho): used with g(rho) in induced dvort
(NULL for unsupported)
- combined(rho): combines g and zeta for perf.
(NULL for unsupported)
- eta(rho): used in particle strength exchange
(NULL for unsupported)
- cl_kernel_name_ext: identifies opencl kernel variant to run.
(fall back to OpenMP)
- 2D and 3D variants
*/
typedef struct {
float(*g_3D)(float rho);
float(*g_2D)(float rho);
float(*zeta_3D)(float rho);
void(*combined_3D)(float rho, float* g, float* zeta);
float(*eta_3D)(float rho);
float(*eta_2D)(float rho);
char cl_kernel_name_ext[32];
} cvtx_VortFunc;

/* cvtx_libary controls */
/* cvtx libary accelerator controls */
CVTX_EXPORT void cvtx_initialise();
CVTX_EXPORT void cvtx_finalise();
CVTX_EXPORT int cvtx_num_accelerators();
Expand All @@ -66,7 +89,7 @@ CVTX_EXPORT int cvtx_accelerator_enabled(int accelerator_id);
CVTX_EXPORT void cvtx_accelerator_enable(int accelerator_id);
CVTX_EXPORT void cvtx_accelerator_disable(int accelerator_id);

/* cvtx_P3D functions */
/* cvtx_P3D 3D vortex particle functions */
CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_vel(
const cvtx_P3D *self,
const bsv_V3f mes_point,
Expand Down Expand Up @@ -142,8 +165,7 @@ CVTX_EXPORT const cvtx_VortFunc cvtx_VortFunc_winckelmans(void);
CVTX_EXPORT const cvtx_VortFunc cvtx_VortFunc_planetary(void);
CVTX_EXPORT const cvtx_VortFunc cvtx_VortFunc_gaussian(void);

/* cvtx_straight vortex filament functions */

/* cvtx_F3D straight vortex filament functions */
CVTX_EXPORT bsv_V3f cvtx_F3D_S2S_vel(
const cvtx_F3D *self,
const bsv_V3f mes_point);
Expand Down Expand Up @@ -184,4 +206,27 @@ CVTX_EXPORT void cvtx_F3D_inf_mtrx(
const int num_mes,
float *result_matrix);

/* cvtx_P2D vortex particle 2D functions */
CVTX_EXPORT bsv_V2f cvtx_P2D_S2S_vel(
const cvtx_P2D *self,
const bsv_V2f mes_point,
const cvtx_VortFunc *kernel,
float regularisation_radius);

CVTX_EXPORT bsv_V2f cvtx_P2D_M2S_vel(
const cvtx_P2D **array_start,
const int num_particles,
const bsv_V2f mes_point,
const cvtx_VortFunc *kernel,
float regularisation_radius);

CVTX_EXPORT void cvtx_P2D_M2M_vel(
const cvtx_P2D **array_start,
const int num_particles,
const bsv_V2f *mes_start,
const int num_mes,
bsv_V2f *result_array,
const cvtx_VortFunc *kernel,
float regularisation_radius);

#endif /* CVTX_LIBCVTX_H */
138 changes: 138 additions & 0 deletions src/P2D.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
#include "libcvtx.h"
/*============================================================================
P2D.c
Vortex particle in 2D with CPU based code.
Copyright(c) 2019 HJA Bird
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files(the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions :
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
============================================================================*/

#include <assert.h>
#include <math.h>
#include <stddef.h>

#ifdef CVTX_USING_OPENCL
# include "ocl_P2D.h"
#endif

/* The induced velocity for a particle excluding the constant
coefficient 1 / 4pi */
inline bsv_V2f P2D_vel_inner(
const cvtx_P2D * self,
const bsv_V2f mes_point,
const cvtx_VortFunc * kernel,
float recip_reg_rad)
{
bsv_V2f rad, ret;
float radd, rho, g;
if (bsv_V2f_isequal(self->coord, mes_point)) {
ret = bsv_V2f_zero();
}
else {
rad = bsv_V2f_minus(mes_point, self->coord);
radd = bsv_V2f_abs(rad);
rho = radd * recip_reg_rad;
g = -kernel->g_2D(rho);
ret.x[0] = rad.x[1] * self->vorticity * g / (radd * radd);
ret.x[1] = -rad.x[0] * self->vorticity * g / (radd * radd);
}
return ret;
}

CVTX_EXPORT bsv_V2f cvtx_P2D_S2S_vel(
const cvtx_P2D * self,
const bsv_V2f mes_point,
const cvtx_VortFunc * kernel,
float regularisation_radius)
{
bsv_V2f ret;
ret = P2D_vel_inner(self, mes_point, kernel,
1.f / fabsf(regularisation_radius));
return bsv_V2f_mult(ret, -1.f / (2.f * acosf(-1.f)));
}

CVTX_EXPORT bsv_V2f cvtx_P2D_M2S_vel(
const cvtx_P2D **array_start,
const int num_particles,
const bsv_V2f mes_point,
const cvtx_VortFunc *kernel,
float regularisation_radius)
{
double rx = 0, ry = 0;
long i;
float recip_reg_rad = 1.f / fabsf(regularisation_radius);
assert(num_particles >= 0);
#pragma omp parallel for reduction(+:rx, ry)
for (i = 0; i < num_particles; ++i) {
bsv_V2f vel = P2D_vel_inner(array_start[i],
mes_point, kernel, recip_reg_rad);
rx += vel.x[0];
ry += vel.x[1];
}
bsv_V2f ret = { (float)rx, (float)ry };
return bsv_V2f_mult(ret, -1.f / (2.f * acosf(-1.f)));
}


static void cpu_brute_force_P2D_M2M_vel(
const cvtx_P2D **array_start,
const int num_particles,
const bsv_V2f *mes_start,
const int num_mes,
bsv_V2f *result_array,
const cvtx_VortFunc *kernel,
float regularisation_radius)
{
long i;
#pragma omp parallel for schedule(static)
for (i = 0; i < num_mes; ++i) {
result_array[i] = cvtx_P2D_M2S_vel(
array_start, num_particles, mes_start[i],
kernel, regularisation_radius);
}
return;
}

CVTX_EXPORT void cvtx_P2D_M2M_vel(
const cvtx_P2D **array_start,
const int num_particles,
const bsv_V2f *mes_start,
const int num_mes,
bsv_V2f *result_array,
const cvtx_VortFunc *kernel,
float regularisation_radius)
{
#ifdef CVTX_USING_OPENCL
if (num_particles < 256
|| num_mes < 256
|| kernel->cl_kernel_name_ext == ""
|| opencl_brute_force_P2D_M2M_vel(
array_start, num_particles, mes_start,
num_mes, result_array, kernel, regularisation_radius) != 0)
#endif
{
cpu_brute_force_P2D_M2M_vel(
array_start, num_particles, mes_start,
num_mes, result_array, kernel, regularisation_radius);
}
return;
}

36 changes: 16 additions & 20 deletions src/P3D.c
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
#include "libcvtx.h"
/*============================================================================
Particle.c
P3D.c
Basic representation of a vortex particle.
Vortex particle in 2D with CPU based code.
Copyright(c) 2018 HJA Bird
Copyright(c) 2019 HJA Bird
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files(the "Software"), to deal
Expand Down Expand Up @@ -35,7 +35,7 @@ SOFTWARE.

/* The induced velocity for a particle excluding the constant
coefficient 1 / 4pi */
inline bsv_V3f particle_ind_vel_inner(
inline bsv_V3f P3D_vel_inner(
const cvtx_P3D * self,
const bsv_V3f mes_point,
const cvtx_VortFunc * kernel,
Expand All @@ -50,7 +50,7 @@ inline bsv_V3f particle_ind_vel_inner(
rad = bsv_V3f_minus(mes_point, self->coord);
radd = bsv_V3f_abs(rad);
rho = radd * recip_reg_rad; /* Assume positive. */
cor = -kernel->g_fn(rho);
cor = -kernel->g_3D(rho);
den = powf(radd, -3);
num = bsv_V3f_cross(rad, self->vorticity);
ret = bsv_V3f_mult(num, cor * den);
Expand All @@ -65,7 +65,7 @@ CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_vel(
float regularisation_radius)
{
bsv_V3f ret;
ret = particle_ind_vel_inner(self, mes_point, kernel,
ret = P3D_vel_inner(self, mes_point, kernel,
1.f/fabsf(regularisation_radius));
return bsv_V3f_mult(ret, 1.f / (4.f * acosf(-1.f)));
}
Expand All @@ -84,7 +84,7 @@ CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_dvort(
rad = bsv_V3f_minus(induced_particle->coord, self->coord);
radd = bsv_V3f_abs(rad);
rho = fabsf(radd / regularisation_radius);
kernel->combined_fn(rho, &g, &f);
kernel->combined_3D(rho, &g, &f);
cross_om = bsv_V3f_cross(induced_particle->vorticity, self->vorticity);
t1 = (float)1. / ((float)4. * (float)acos(-1) * powf(regularisation_radius, 3));
t21n = bsv_V3f_mult(cross_om, g);
Expand All @@ -100,10 +100,6 @@ CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_dvort(
return ret;
}

float sphere_volume(float radius){
return 4 * (float)acos(-1) * radius * radius * radius / (float) 3.;
}

CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_visc_dvort(
const cvtx_P3D * self,
const cvtx_P3D * induced_particle,
Expand All @@ -113,7 +109,7 @@ CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_visc_dvort(
{
bsv_V3f ret, rad, t211, t212, t21, t2;
float radd, rho, t1, t22;
assert(kernel->eta_fn != NULL && "Used vortex regularisation"
assert(kernel->eta_3D != NULL && "Used vortex regularisation"
"that did have a defined eta function");
if(bsv_V3f_isequal(self->coord, induced_particle->coord)){
ret = bsv_V3f_zero();
Expand All @@ -128,7 +124,7 @@ CVTX_EXPORT bsv_V3f cvtx_P3D_S2S_visc_dvort(
t212 = bsv_V3f_mult(induced_particle->vorticity,
-1 * self->volume);
t21 = bsv_V3f_plus(t211, t212);
t22 = kernel->eta_fn(rho);
t22 = kernel->eta_3D(rho);
t2 = bsv_V3f_mult(t21, t22);
ret = bsv_V3f_mult(t2, t1);
}
Expand All @@ -148,7 +144,7 @@ CVTX_EXPORT bsv_V3f cvtx_P3D_M2S_vel(
assert(num_particles >= 0);
#pragma omp parallel for reduction(+:rx, ry, rz)
for (i = 0; i < num_particles; ++i) {
bsv_V3f vel = particle_ind_vel_inner(array_start[i],
bsv_V3f vel = P3D_vel_inner(array_start[i],
mes_point, kernel, recip_reg_rad);
rx += vel.x[0];
ry += vel.x[1];
Expand Down Expand Up @@ -203,7 +199,7 @@ CVTX_EXPORT bsv_V3f cvtx_P3D_M2S_visc_dvort(
return ret;
}

static void cpu_brute_force_ParticleArr_Arr_ind_vel(
static void cpu_brute_force_P3D_M2M_vel(
const cvtx_P3D **array_start,
const int num_particles,
const bsv_V3f *mes_start,
Expand Down Expand Up @@ -240,14 +236,14 @@ CVTX_EXPORT void cvtx_P3D_M2M_vel(
num_mes, result_array, kernel, regularisation_radius) != 0)
#endif
{
cpu_brute_force_ParticleArr_Arr_ind_vel(
cpu_brute_force_P3D_M2M_vel(
array_start, num_particles, mes_start,
num_mes, result_array, kernel, regularisation_radius);
}
return;
}

void cpu_brute_force_ParticleArr_Arr_ind_dvort(
void cpu_brute_force_P3D_M2M_dvort(
const cvtx_P3D **array_start,
const int num_particles,
const cvtx_P3D **induced_start,
Expand Down Expand Up @@ -284,14 +280,14 @@ CVTX_EXPORT void cvtx_P3D_M2M_dvort(
num_induced, result_array, kernel, regularisation_radius) != 0)
#endif
{
cpu_brute_force_ParticleArr_Arr_ind_dvort(
cpu_brute_force_P3D_M2M_dvort(
array_start, num_particles, induced_start,
num_induced, result_array, kernel, regularisation_radius);
}
return;
}

void cpu_brute_force_ParticleArr_Arr_visc_ind_dvort(
void cpu_brute_force_P3D_M2M_visc_dvort(
const cvtx_P3D **array_start,
const int num_particles,
const cvtx_P3D **induced_start,
Expand Down Expand Up @@ -330,7 +326,7 @@ CVTX_EXPORT void cvtx_P3D_M2M_visc_dvort(
num_induced, result_array, kernel, regularisation_radius, kinematic_visc) != 0)
#endif
{
cpu_brute_force_ParticleArr_Arr_visc_ind_dvort(
cpu_brute_force_P3D_M2M_visc_dvort(
array_start, num_particles, induced_start,
num_induced, result_array, kernel, regularisation_radius, kinematic_visc);
}
Expand Down
Loading

0 comments on commit 2017561

Please sign in to comment.