Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
TurTLE
TurTLE
Commits
6814dd0c
Commit
6814dd0c
authored
Dec 26, 2015
by
Cristian Lalescu
Browse files
Merge branch 'feature/rFFTW_optimizations' into develop
parents
83bba017
fffb6965
Changes
14
Hide whitespace changes
Inline
Side-by-side
.gitignore
View file @
6814dd0c
...
...
@@ -6,3 +6,7 @@ meta/.ipynb_checkpoints/*
dist
build
bfps.egg-info
MANIFEST.in
bfps/install_info.pickle
bfps/lib*.a
obj
bfps/NavierStokes.py
View file @
6814dd0c
...
...
@@ -408,8 +408,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
update_fields
+=
'fs->low_pass_Fourier(fs->cvelocity, 3, {0});
\n
'
.
format
(
kcut
)
update_fields
+=
(
'fs->ift_velocity();
\n
'
+
'vel_{0}->read_rFFTW(fs->rvelocity);
\n
'
+
'fs->compute_Lagrangian_acceleration(acc_{0}->temp);
\n
'
+
'acc_{0}->read_rFFTW(acc_{0}->temp);
\n
'
).
format
(
name
)
'fs->compute_Lagrangian_acceleration(acc_{0}->field);
\n
'
).
format
(
name
)
self
.
fluid_start
+=
update_fields
self
.
fluid_loop
+=
update_fields
return
None
...
...
@@ -450,7 +449,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
double *velocity = new double[ps{0}->array_size];
ps{0}->sample_vec_field(vel_{1}, velocity);
ps{0}->sample_vec_field(acc_{1}, acceleration);
if (ps{0}->
fs->rd->
myrank == 0)
if (ps{0}->myrank == 0)
{{
//VELOCITY begin
std::string temp_string = (std::string("/particles/") +
...
...
@@ -512,7 +511,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
neighbours
,
self
.
particle_species
)
self
.
particle_start
+=
(
'ps{0} = new rFFTW_particles<VELOCITY_TRACER, {1}, {2}>(
\n
'
+
'fname,
fs,
vel_{3},
\n
'
+
'fname, vel_{3},
\n
'
+
'nparticles,
\n
'
+
'niter_part, tracers{0}_integration_steps);
\n
'
).
format
(
self
.
particle_species
,
self
.
C_dtype
,
neighbours
,
fields_name
)
...
...
@@ -536,10 +535,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self
.
particle_loop
+=
'ps{0}->synchronize();
\n
'
.
format
(
self
.
particle_species
)
elif
particle_class
==
'rFFTW_particles'
:
self
.
particle_loop
+=
'ps{0}->step();
\n
'
.
format
(
self
.
particle_species
)
self
.
particle_loop
+=
((
'if (ps{0}->iteration % niter_part == 0)
\n
'
+
'{{
\n
'
+
'ps{0}->write(stat_file, false);
\n
'
).
format
(
self
.
particle_species
)
+
output_vel_acc
+
'}
\n
'
)
self
.
particle_stat_src
+=
(
(
'if (ps{0}->iteration % niter_part == 0)
\n
'
+
'{{
\n
'
+
'ps{0}->write(stat_file, false);
\n
'
).
format
(
self
.
particle_species
)
+
output_vel_acc
+
'}
\n
'
)
self
.
particle_species
+=
1
return
None
def
get_data_file
(
self
):
...
...
bfps/cpp/rFFTW_interpolator.cpp
View file @
6814dd0c
...
...
@@ -26,6 +26,7 @@
#define NDEBUG
#include
<cmath>
#include
"rFFTW_interpolator.hpp"
template
<
class
rnumber
,
int
interp_neighbours
>
...
...
@@ -37,58 +38,104 @@ rFFTW_interpolator<rnumber, interp_neighbours>::rFFTW_interpolator(
this
->
field_size
=
2
*
fs
->
cd
->
local_size
;
this
->
compute_beta
=
BETA_POLYS
;
if
(
sizeof
(
rnumber
)
==
4
)
{
this
->
f0
=
(
rnumber
*
)((
void
*
)
fftwf_alloc_real
(
this
->
field_size
));
this
->
f1
=
(
rnumber
*
)((
void
*
)
fftwf_alloc_real
(
this
->
field_size
));
}
this
->
field
=
(
rnumber
*
)((
void
*
)
fftwf_alloc_real
(
this
->
field_size
));
else
if
(
sizeof
(
rnumber
)
==
8
)
{
this
->
f0
=
(
rnumber
*
)((
void
*
)
fftw_alloc_real
(
this
->
field_size
));
this
->
f1
=
(
rnumber
*
)((
void
*
)
fftw_alloc_real
(
this
->
field_size
));
}
this
->
temp
=
this
->
f1
;
this
->
field
=
(
rnumber
*
)((
void
*
)
fftw_alloc_real
(
this
->
field_size
));
// compute dx, dy, dz;
this
->
dx
=
4
*
acos
(
0
)
/
(
fs
->
dkx
*
this
->
descriptor
->
sizes
[
2
]);
this
->
dy
=
4
*
acos
(
0
)
/
(
fs
->
dky
*
this
->
descriptor
->
sizes
[
1
]);
this
->
dz
=
4
*
acos
(
0
)
/
(
fs
->
dkz
*
this
->
descriptor
->
sizes
[
0
]);
// generate compute array
this
->
compute
=
new
bool
[
this
->
descriptor
->
sizes
[
0
]];
std
::
fill_n
(
this
->
compute
,
this
->
descriptor
->
sizes
[
0
],
false
);
for
(
int
iz
=
this
->
descriptor
->
starts
[
0
]
-
interp_neighbours
-
1
;
iz
<=
this
->
descriptor
->
starts
[
0
]
+
this
->
descriptor
->
subsizes
[
0
]
+
interp_neighbours
;
iz
++
)
this
->
compute
[((
iz
+
this
->
descriptor
->
sizes
[
0
])
%
this
->
descriptor
->
sizes
[
0
])]
=
true
;
}
template
<
class
rnumber
,
int
interp_neighbours
>
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>::~
rFFTW_interpolator
()
{
if
(
sizeof
(
rnumber
)
==
4
)
{
fftwf_free
((
float
*
)((
void
*
)
this
->
f0
));
fftwf_free
((
float
*
)((
void
*
)
this
->
f1
));
}
fftwf_free
((
float
*
)((
void
*
)
this
->
field
));
else
if
(
sizeof
(
rnumber
)
==
8
)
{
fftw_free
((
double
*
)((
void
*
)
this
->
f0
));
fftw_free
((
double
*
)((
void
*
)
this
->
f1
));
}
fftw_free
((
double
*
)((
void
*
)
this
->
field
));
delete
[]
this
->
compute
;
}
template
<
class
rnumber
,
int
interp_neighbours
>
int
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>::
read_rFFTW
(
void
*
void_src
)
{
/* first, roll fields */
rnumber
*
tmp
=
this
->
f0
;
this
->
f0
=
this
->
f1
;
this
->
f1
=
tmp
;
this
->
temp
=
this
->
f0
;
/* now do regular things */
rnumber
*
src
=
(
rnumber
*
)
void_src
;
rnumber
*
dst
=
this
->
f1
;
/* do big copy of middle stuff */
std
::
copy
(
src
,
src
+
this
->
field_size
,
dst
);
this
->
field
);
return
EXIT_SUCCESS
;
}
template
<
class
rnumber
,
int
interp_neighbours
>
void
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>::
operator
()(
double
t
,
void
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>::
get_grid_coordinates
(
const
int
nparticles
,
const
int
pdimension
,
const
double
*
x
,
int
*
xg
,
double
*
xx
,
double
*
xx
)
{
static
double
grid_size
[]
=
{
this
->
dx
,
this
->
dy
,
this
->
dz
};
double
tval
;
std
::
fill_n
(
xg
,
nparticles
*
3
,
0
);
std
::
fill_n
(
xx
,
nparticles
*
3
,
0.0
);
for
(
int
p
=
0
;
p
<
nparticles
;
p
++
)
{
for
(
int
c
=
0
;
c
<
3
;
c
++
)
{
tval
=
floor
(
x
[
p
*
pdimension
+
c
]
/
grid_size
[
c
]);
xg
[
p
*
3
+
c
]
=
MOD
(
int
(
tval
),
this
->
descriptor
->
sizes
[
2
-
c
]);
xx
[
p
*
3
+
c
]
=
(
x
[
p
*
pdimension
+
c
]
-
tval
*
grid_size
[
c
])
/
grid_size
[
c
];
}
}
}
template
<
class
rnumber
,
int
interp_neighbours
>
void
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>::
sample
(
const
int
nparticles
,
const
int
pdimension
,
const
double
*
__restrict__
x
,
double
*
__restrict__
y
,
const
int
*
deriv
)
{
/* get grid coordinates */
int
*
xg
=
new
int
[
3
*
nparticles
];
double
*
xx
=
new
double
[
3
*
nparticles
];
double
*
yy
=
new
double
[
3
*
nparticles
];
std
::
fill_n
(
yy
,
3
*
nparticles
,
0.0
);
this
->
get_grid_coordinates
(
nparticles
,
pdimension
,
x
,
xg
,
xx
);
/* perform interpolation */
for
(
int
p
=
0
;
p
<
nparticles
;
p
++
)
if
(
this
->
compute
[
xg
[
p
*
3
+
2
]])
this
->
operator
()(
xg
+
p
*
3
,
xx
+
p
*
3
,
yy
+
p
*
3
,
deriv
);
MPI_Allreduce
(
yy
,
y
,
3
*
nparticles
,
MPI_DOUBLE
,
MPI_SUM
,
this
->
descriptor
->
comm
);
delete
[]
yy
;
delete
[]
xg
;
delete
[]
xx
;
}
template
<
class
rnumber
,
int
interp_neighbours
>
void
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>::
operator
()(
const
int
*
xg
,
const
double
*
xx
,
double
*
dest
,
int
*
deriv
)
const
int
*
deriv
)
{
double
bx
[
interp_neighbours
*
2
+
2
],
by
[
interp_neighbours
*
2
+
2
],
bz
[
interp_neighbours
*
2
+
2
];
if
(
deriv
==
NULL
)
...
...
@@ -105,13 +152,11 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
}
std
::
fill_n
(
dest
,
3
,
0
);
ptrdiff_t
bigiz
,
bigiy
,
bigix
;
double
tval
[
3
];
for
(
int
iz
=
-
interp_neighbours
;
iz
<=
interp_neighbours
+
1
;
iz
++
)
{
bigiz
=
ptrdiff_t
(((
xg
[
2
]
+
iz
)
+
this
->
descriptor
->
sizes
[
0
])
%
this
->
descriptor
->
sizes
[
0
]);
if
(
this
->
descriptor
->
myrank
==
this
->
descriptor
->
rank
[
bigiz
])
{
std
::
fill_n
(
tval
,
3
,
0
);
for
(
int
iy
=
-
interp_neighbours
;
iy
<=
interp_neighbours
+
1
;
iy
++
)
{
bigiy
=
ptrdiff_t
(
MOD
(
xg
[
1
]
+
iy
,
this
->
descriptor
->
sizes
[
1
]));
...
...
@@ -122,17 +167,11 @@ void rFFTW_interpolator<rnumber, interp_neighbours>::operator()(
bigiy
)
*
(
this
->
descriptor
->
sizes
[
2
]
+
2
)
+
bigix
)
*
3
;
for
(
int
c
=
0
;
c
<
3
;
c
++
)
{
dest
[
c
]
+=
(
this
->
f0
[
tindex
+
c
]
*
(
1
-
t
)
+
t
*
this
->
f1
[
tindex
+
c
])
*
(
bz
[
iz
+
interp_neighbours
]
*
by
[
iy
+
interp_neighbours
]
*
bx
[
ix
+
interp_neighbours
]);
tval
[
c
]
+=
(
this
->
f0
[
tindex
+
c
]
*
(
1
-
t
)
+
t
*
this
->
f1
[
tindex
+
c
])
*
(
bz
[
iz
+
interp_neighbours
]
*
by
[
iy
+
interp_neighbours
]
*
bx
[
ix
+
interp_neighbours
]);
}
dest
[
c
]
+=
this
->
field
[
tindex
+
c
]
*
(
bz
[
iz
+
interp_neighbours
]
*
by
[
iy
+
interp_neighbours
]
*
bx
[
ix
+
interp_neighbours
]);
}
}
DEBUG_MSG
(
"%ld %d %d %g %g %g
\n
"
,
bigiz
,
xg
[
1
],
xg
[
0
],
tval
[
0
],
tval
[
1
],
tval
[
2
]);
}
}
}
...
...
bfps/cpp/rFFTW_interpolator.hpp
View file @
6814dd0c
...
...
@@ -44,17 +44,52 @@ template <class rnumber, int interp_neighbours>
class
rFFTW_interpolator
{
public:
/* size of field to interpolate */
ptrdiff_t
field_size
;
/* pointer to polynomial function */
base_polynomial_values
compute_beta
;
/* descriptor of field to interpolate */
field_descriptor
<
rnumber
>
*
descriptor
;
rnumber
*
f0
,
*
f1
,
*
temp
;
/* pointers to fields that are to be interpolated
* */
rnumber
*
field
;
/* physical parameters of field */
double
dx
,
dy
,
dz
;
/* compute[iz] is true if .
* local_zstart - neighbours <= iz <= local_zend + 1 + neighbours
* */
bool
*
compute
;
rFFTW_interpolator
(
fluid_solver_base
<
rnumber
>
*
FSOLVER
,
base_polynomial_values
BETA_POLYS
);
~
rFFTW_interpolator
();
void
operator
()(
double
t
,
int
*
__restrict__
xg
,
double
*
__restrict__
xx
,
double
*
__restrict__
dest
,
int
*
deriv
=
NULL
);
/* map real locations to grid coordinates */
void
get_grid_coordinates
(
const
int
nparticles
,
const
int
pdimension
,
const
double
*
__restrict__
x
,
int
*
__restrict__
xg
,
double
*
__restrict__
xx
);
/* interpolate field at an array of locations */
void
sample
(
const
int
nparticles
,
const
int
pdimension
,
const
double
*
__restrict__
x
,
double
*
__restrict__
y
,
const
int
*
deriv
=
NULL
);
/* interpolate 1 point */
void
operator
()(
const
int
*
__restrict__
xg
,
const
double
*
__restrict__
xx
,
double
*
__restrict__
dest
,
const
int
*
deriv
=
NULL
);
int
read_rFFTW
(
void
*
src
);
};
...
...
bfps/cpp/rFFTW_particles.cpp
View file @
6814dd0c
...
...
@@ -24,7 +24,7 @@
//
#define NDEBUG
#define NDEBUG
#include
<cmath>
#include
<cassert>
...
...
@@ -42,7 +42,6 @@ extern int myrank, nprocs;
template
<
int
particle_type
,
class
rnumber
,
int
interp_neighbours
>
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
rFFTW_particles
(
const
char
*
NAME
,
fluid_solver_base
<
rnumber
>
*
FSOLVER
,
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>
*
FIELD
,
const
int
NPARTICLES
,
const
int
TRAJ_SKIP
,
...
...
@@ -60,7 +59,6 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
assert
((
INTEGRATION_STEPS
<=
6
)
&&
(
INTEGRATION_STEPS
>=
1
));
strncpy
(
this
->
name
,
NAME
,
256
);
this
->
fs
=
FSOLVER
;
this
->
nparticles
=
NPARTICLES
;
this
->
vel
=
FIELD
;
this
->
integration_steps
=
INTEGRATION_STEPS
;
...
...
@@ -76,41 +74,6 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::rFFTW_particles(
this
->
rhs
[
i
]
=
new
double
[
this
->
array_size
];
std
::
fill_n
(
this
->
rhs
[
i
],
this
->
array_size
,
0.0
);
}
// compute dx, dy, dz;
this
->
dx
=
4
*
acos
(
0
)
/
(
this
->
fs
->
dkx
*
this
->
fs
->
rd
->
sizes
[
2
]);
this
->
dy
=
4
*
acos
(
0
)
/
(
this
->
fs
->
dky
*
this
->
fs
->
rd
->
sizes
[
1
]);
this
->
dz
=
4
*
acos
(
0
)
/
(
this
->
fs
->
dkz
*
this
->
fs
->
rd
->
sizes
[
0
]);
// compute lower and upper bounds
this
->
lbound
=
new
double
[
nprocs
];
this
->
ubound
=
new
double
[
nprocs
];
double
*
tbound
=
new
double
[
nprocs
];
std
::
fill_n
(
tbound
,
nprocs
,
0.0
);
tbound
[
this
->
myrank
]
=
this
->
fs
->
rd
->
starts
[
0
]
*
this
->
dz
;
MPI_Allreduce
(
tbound
,
this
->
lbound
,
nprocs
,
MPI_DOUBLE
,
MPI_SUM
,
this
->
comm
);
std
::
fill_n
(
tbound
,
nprocs
,
0.0
);
tbound
[
this
->
myrank
]
=
(
this
->
fs
->
rd
->
starts
[
0
]
+
this
->
fs
->
rd
->
subsizes
[
0
])
*
this
->
dz
;
MPI_Allreduce
(
tbound
,
this
->
ubound
,
nprocs
,
MPI_DOUBLE
,
MPI_SUM
,
this
->
comm
);
delete
[]
tbound
;
//for (int r = 0; r<nprocs; r++)
// DEBUG_MSG(
// "lbound[%d] = %lg, ubound[%d] = %lg\n",
// r, this->lbound[r],
// r, this->ubound[r]
// );
}
template
<
int
particle_type
,
class
rnumber
,
int
interp_neighbours
>
...
...
@@ -121,59 +84,19 @@ rFFTW_particles<particle_type, rnumber, interp_neighbours>::~rFFTW_particles()
{
delete
[]
this
->
rhs
[
i
];
}
delete
[]
this
->
lbound
;
delete
[]
this
->
ubound
;
}
template
<
int
particle_type
,
class
rnumber
,
int
interp_neighbours
>
void
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
sample_vec_field
(
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>
*
vec
,
double
t
,
double
*
x
,
double
*
y
,
int
*
deriv
)
{
/* get grid coordinates */
int
*
xg
=
new
int
[
3
*
this
->
nparticles
];
double
*
xx
=
new
double
[
3
*
this
->
nparticles
];
double
*
yy
=
new
double
[
3
*
this
->
nparticles
];
std
::
fill_n
(
yy
,
3
*
this
->
nparticles
,
0.0
);
this
->
get_grid_coordinates
(
x
,
xg
,
xx
);
/* perform interpolation */
for
(
int
p
=
0
;
p
<
this
->
nparticles
;
p
++
)
(
*
vec
)(
t
,
xg
+
p
*
3
,
xx
+
p
*
3
,
yy
+
p
*
3
,
deriv
);
MPI_Allreduce
(
yy
,
y
,
3
*
this
->
nparticles
,
MPI_DOUBLE
,
MPI_SUM
,
this
->
comm
);
delete
[]
yy
;
delete
[]
xg
;
delete
[]
xx
;
}
template
<
int
particle_type
,
class
rnumber
,
int
interp_neighbours
>
void
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
get_rhs
(
double
t
,
double
*
x
,
double
*
y
)
void
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
get_rhs
(
double
*
x
,
double
*
y
)
{
switch
(
particle_type
)
{
case
VELOCITY_TRACER
:
this
->
sample
_vec_field
(
this
->
vel
,
t
,
x
,
y
);
this
->
vel
->
sample
(
this
->
nparticles
,
this
->
ncomponents
,
x
,
y
);
break
;
}
}
template
<
int
particle_type
,
class
rnumber
,
int
interp_neighbours
>
int
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
get_rank
(
double
z
)
{
int
tmp
=
this
->
fs
->
rd
->
rank
[
MOD
(
int
(
floor
(
z
/
this
->
dz
)),
this
->
fs
->
rd
->
sizes
[
0
])];
assert
(
tmp
>=
0
&&
tmp
<
this
->
fs
->
rd
->
nprocs
);
return
tmp
;
}
template
<
int
particle_type
,
class
rnumber
,
int
interp_neighbours
>
void
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
roll_rhs
()
{
...
...
@@ -189,7 +112,7 @@ template <int particle_type, class rnumber, int interp_neighbours>
void
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
AdamsBashforth
(
int
nsteps
)
{
ptrdiff_t
ii
;
this
->
get_rhs
(
0
,
this
->
state
,
this
->
rhs
[
0
]);
this
->
get_rhs
(
this
->
state
,
this
->
rhs
[
0
]);
switch
(
nsteps
)
{
case
1
:
...
...
@@ -270,29 +193,6 @@ void rFFTW_particles<particle_type, rnumber, interp_neighbours>::step()
}
template
<
int
particle_type
,
class
rnumber
,
int
interp_neighbours
>
void
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
get_grid_coordinates
(
double
*
x
,
int
*
xg
,
double
*
xx
)
{
static
double
grid_size
[]
=
{
this
->
dx
,
this
->
dy
,
this
->
dz
};
double
tval
;
std
::
fill_n
(
xg
,
this
->
nparticles
*
3
,
0
);
std
::
fill_n
(
xx
,
this
->
nparticles
*
3
,
0.0
);
for
(
int
p
=
0
;
p
<
this
->
nparticles
;
p
++
)
{
for
(
int
c
=
0
;
c
<
3
;
c
++
)
{
tval
=
floor
(
x
[
p
*
this
->
ncomponents
+
c
]
/
grid_size
[
c
]);
xg
[
p
*
3
+
c
]
=
MOD
(
int
(
tval
),
this
->
fs
->
rd
->
sizes
[
2
-
c
]);
xx
[
p
*
3
+
c
]
=
(
x
[
p
*
this
->
ncomponents
+
c
]
-
tval
*
grid_size
[
c
])
/
grid_size
[
c
];
}
/*xg[p*3+2] -= this->fs->rd->starts[0];
if (this->myrank == this->fs->rd->rank[0] &&
xg[p*3+2] > this->fs->rd->subsizes[0])
xg[p*3+2] -= this->fs->rd->sizes[0];*/
}
}
template
<
int
particle_type
,
class
rnumber
,
int
interp_neighbours
>
void
rFFTW_particles
<
particle_type
,
rnumber
,
interp_neighbours
>::
read
(
hid_t
data_file_id
)
{
...
...
bfps/cpp/rFFTW_particles.hpp
View file @
6814dd0c
...
...
@@ -43,20 +43,8 @@ class rFFTW_particles
public:
int
myrank
,
nprocs
;
MPI_Comm
comm
;
fluid_solver_base
<
rnumber
>
*
fs
;
rnumber
*
data
;
/* watching is an array of shape [nparticles], with
* watching[p] being true if particle p is in the domain of myrank
* or in the buffer regions.
* only used if multistep is false.
* */
bool
*
watching
;
/* computing is an array of shape [nparticles], with
* computing[p] being the rank that is currently working on particle p
* */
int
*
computing
;
/* state will generally hold all the information about the particles.
* in the beginning, we will only need to solve 3D ODEs, but I figured
* a general ncomponents is better, since we may change our minds.
...
...
@@ -68,8 +56,6 @@ class rFFTW_particles
int
array_size
;
int
integration_steps
;
int
traj_skip
;
double
*
lbound
;
double
*
ubound
;
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>
*
vel
;
/* simulation parameters */
...
...
@@ -77,23 +63,15 @@ class rFFTW_particles
int
iteration
;
double
dt
;
/* physical parameters of field */
double
dx
,
dy
,
dz
;
/* methods */
/* constructor and destructor.
* allocate and deallocate:
* this->state
* this->rhs
* this->lbound
* this->ubound
* this->computing
* this->watching
* */
rFFTW_particles
(
const
char
*
NAME
,
fluid_solver_base
<
rnumber
>
*
FSOLVER
,
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>
*
FIELD
,
const
int
NPARTICLES
,
const
int
TRAJ_SKIP
,
...
...
@@ -101,19 +79,10 @@ class rFFTW_particles
~
rFFTW_particles
();
void
get_rhs
(
double
*
__restrict__
x
,
double
*
__restrict__
rhs
);
void
get_rhs
(
double
t
,
double
*
__restrict__
x
,
double
*
__restrict__
rhs
);
int
get_rank
(
double
z
);
// get rank for given value of z
void
get_grid_coordinates
(
double
*
__restrict__
x
,
int
*
__restrict__
xg
,
double
*
__restrict__
xx
);
void
sample_vec_field
(
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>
*
vec
,
double
t
,
double
*
__restrict__
x
,
double
*
__restrict__
y
,
int
*
deriv
=
NULL
);
inline
void
sample_vec_field
(
rFFTW_interpolator
<
rnumber
,
interp_neighbours
>
*
field
,
double
*
vec_values
)
{
this
->
sample
_vec_field
(
field
,
1.0
,
this
->
state
,
vec_values
,
NULL
);
field
->
sample
(
this
->
nparticles
,
this
->
ncomponents
,
this
->
state
,
vec_values
,
NULL
);
}
/* input/output */
...
...
bfps/fluid_base.py
View file @
6814dd0c
...
...
@@ -90,13 +90,14 @@ class fluid_particle_base(bfps.code):
self
.
fluid_loop
=
''
self
.
fluid_end
=
''
self
.
fluid_output
=
''
self
.
stat_src
=
''
self
.
particle_includes
=
''
self
.
particle_variables
=
''
self
.
particle_definitions
=
''
self
.
particle_start
=
''
self
.
particle_loop
=
''
self
.
particle_end
=
''
self
.
stat_src
=
''
self
.
particle_
stat_src
=
''
self
.
file_datasets_grow
=
''
return
None
def
finalize_code
(
self
):
...
...
@@ -144,6 +145,7 @@ class fluid_particle_base(bfps.code):
'return file_problems;
\n
'
'}
\n
'
)
self
.
definitions
+=
'void do_stats()
\n
{
\n
'
+
self
.
stat_src
+
'}
\n
'
self
.
definitions
+=
'void do_particle_stats()
\n
{
\n
'
+
self
.
particle_stat_src
+
'}
\n
'
# take care of wisdom
if
self
.
use_fftw_wisdom
:
if
self
.
dtype
==
np
.
float32
:
...
...
@@ -191,6 +193,7 @@ class fluid_particle_base(bfps.code):
return EXIT_SUCCESS;
}
do_stats();
do_particle_stats();
//endcpp
"""
output_time_difference
=
(
'time1 = clock();
\n
'
+
...
...
@@ -202,16 +205,21 @@ class fluid_particle_base(bfps.code):
'<< iteration << " took " '
+
'<< time_difference/nprocs << " seconds" << std::endl;
\n
'
+
'time0 = time1;
\n
'
)
self
.
main
+=
'for (int max_iter = iteration+niter_todo; iteration < max_iter; iteration++)
\n
{
\n
'
self
.
main
+=
'for (int max_iter = iteration+niter_todo; iteration < max_iter; iteration++)
\n
'
self
.
main
+=
'{
\n
'
self
.
main
+=
output_time_difference
self
.
main
+=
self
.
fluid_loop
if
self
.
particle_species
>
0
:
self
.
main
+=
self
.
particle_loop
self
.
main
+=
'if (iteration % niter_stat == 0) do_stats();
\n
}
\n
'
self
.
main
+=
self
.
fluid_loop
self
.
main
+=
'if (iteration % niter_stat == 0) do_stats();
\n
'
if
self
.
particle_species
>
0
:
self
.
main
+=
'if (iteration % niter_part == 0) do_particle_stats();
\n
'
self
.
main
+=
'}
\n
'
self
.
main
+=
output_time_difference
self
.
main
+=
'do_stats();
\n
'
self
.
main
+=
'do_particle_stats();
\n
'
if
self
.
particle_species
>
0
:
self
.
main
+=
self
.
particle_end
self
.
main
+=
'do_stats();
\n
'
self
.
main
+=
self
.
fluid_end
return
None
def
read_rfield
(
...
...
done.txt
View file @
6814dd0c
x 2015-12-04 make code py3 compatible @python3
x 2015-12-23 decide on versioning system +merge0
x 2015-12-24 move get_grid coords to interpolator @optimization +v1.0
x 2015-12-25 get rid of temporal interpolation @optimization +v1.0
x 2015-12-26 call interpolation only when needed @optimization +v1.0
tests/
test_
base.py
→
tests/base.py
View file @
6814dd0c
...
...
@@ -25,9 +25,13 @@
import
os
import
sys
import
subprocess
import
pickle
import
numpy
as
np
import
matplotlib.pyplot
as
plt
import
bfps
from
bfps
import
fluid_resize
...
...
@@ -211,6 +215,93 @@ def acceleration_test(c, m = 3, species = 0):