We present a highly-parallel implementation of the Langevin simulation method for modeling ferrofluids
on Graphical Processor Units (GPU). Our method is based on the Barnes-Hut algorithm. As a benchmark
we use the straightforward 'all-to-all interaction' algorithm. The obtained results are in good agreement
with known theoretical model. With the proposed method we were able to follow the evolution of a
system of one million interacting particles over long time-scales, the task hitherto is out of reach with the
standard, CPU-based numerical schemes.
When you are citing the document, use the following link http://essuir.sumdu.edu.ua/handle/123456789/3534

Denisov, S.

Hanggi, P.

Liutyi, Taras Volodymyrovych

Poliakov, Oleksandr Yuriiovych

Лютий, Тарас Володимирович

Лютый, Тарас Владимирович

Поляков, Александр Юрьевич

Поляков, Олександр Юрійович

SocioEconomic Challenges Journal (Sumy State University)

PROCEEDINGS OF THE INTERNATIONAL CONFERENCE NANOMATERIALS: APPLICATIONS AND PROPERTIES 
Vol. 1 No 1, 04MFPN05(3pp) (2012) 
 
 
2304-1862/2012/1(1)04MFPN05(3) 04MFPN05-1  2012 Sumy State University 
Technique of the Fast Ferrofluids Simulation 
 
A.Yu. Polyakov1,2,*, T.V. Lyutyy1,†, S. Denisov2, P. Hänggi2 
 
1 Sumy State University, 2, Rymsky Korsakov Str., 40007 Sumy, Ukraine 
2 Institut für Physik, Universität Augsburg, Universitätsstraβe 1, D-86135 Augsburg, Germany 
 
(Received 10 July 2012; published online 22 August 2012) 
 
We present a highly-parallel implementation of the Langevin simulation method for modeling ferroflu-
ids on Graphical Processor Units (GPU). Our method is based on the Barnes-Hut algorithm. As a bench-
mark we use the straightforward 'all-to-all interaction' algorithm. The obtained results are in good agree-
ment with known theoretical model. With the proposed method we were able to follow the evolution of a 
system of one million interacting particles over long time-scales, the task hitherto is out of reach with the 
standard, CPU-based numerical schemes. 
 
Keywords: Ferrofluids, Molecular dynamics, N-body problem, Graphics processing unit, CUDA, Barnes-
Hut algorithm. 
 
 PACS numbers: 47.65.Cb, 02.70.Ns 
 
 
                                                          
* alexander.p1987@gmail.com 
† lyutyy@oeph.sumdu.edu.ua 
1. INTRODUCTION 
 
A system of N interacting particles is one of the 
most common set-ups in physics. Since analytical solu-
tions typically do not exist when the number of inter-
acted particles is larger then two, N > 2, the various 
numerical methods are used in order to get an insight 
in the dynamics of the corresponding system. The in-
teraction involves all the pairs of particles, and there-
fore the computational time and efforts scale like N2. 
This effect certainly constitutes a problem, and, in 
order to run big systems, researchers relied on big 
computers. Practically that means the use of large 
computational clusters consisting of many Central 
Processing Units (CPUs). Such clusters are physically 
big, expensive and consume a lot of electricity. Over the 
last decade, an alternative, Graphic Processing Units 
(GPUs), become more and more popular. GPUs, initial-
ly designed to serve as the data-pipelines for graphical 
information, posses a high level of parallelism, while 
being cheap (as compared to CPUs) and small. The 
high computation capability of GPUs is owing to their 
architecture, which represents a set of hundreds of 
computational cores (note that the best modern CPUs 
have typically up to 16 cores), which can perform com-
putations in parallel. Nowadays, the scientific GPU-
computing is used in many areas of computational 
physics, thanks to the Compute United Device Archi-
tecture (CUDA) developed by NVIDIA Corporation, 
which step simplified the development of programs for 
GPUs. 
 
2. MODEL OF THE DISPERSED FERROFLUID 
 
Theory of the continuous bulk ferrofluid is well-
developed within the concept of ferrohydrodynamics [1]. 
However, recently the aspects of the dynamics of confined 
ferrofluids attracted a lot of attention because of the 
promising potential of ferrofluids for biotechnologies and 
medicine [2]. In this context, the disperse structure of a 
ferrofluidic media certainly deserves an attention. The 
dynamics of a ferromagnetic particle, suspended in liq-
uid, can be described by a Langevin equation [3] that 
accounts for the long-range dipole interactions with other 
particles, caused by the particle magnetization, and for 
short-range steric repulsions, caused by the anti-
agglomeration coating. Here we use the simplest model 
and consider an ensemble of identical spherical particles 
of radius R, submerged into the liquid of viscosity . We 
suppose also that the hydrodynamic radius is equal to R, 
and every particle has a homogenous structure with 
respect to its inertia and magnetic properties. Therefore, 
the dynamics of i-th particle is determined by the follow-
ing equations 
 
 
22
2 2
2
,
5
i i ch
ix iy iy ix r r
Td d
u h u h
d d VR D
 (1) 
 
2
2
2
2
2
( )sin
5
( )cos ,
i
iy iz iz iy k
i ch
iz ix ix iz k r r
d
u h u h
d
Td
u h u h
d VR D
 (2) 
 
 
22
2 2
,i idip chsr t di i
Td d
d d VR D
ρ ρ
f f ξ  (3) 
 
where i and  i are the azimuthal and polar angles of i-th 
particle;   t/Tch is the dimensionless time (Tch  
 R/ (3D/4 0)0.5;  is the particle material magnetization, 
0  4 10 – 7H/m is the magnetic constant, D is the 
material particle density); ui  mi/m is the dimension-
less magnetic moment (m  4 R3 /3  V , V is the par-
ticle volume); Γr  8 RTch/DV is the rotation friction 
coefficient; i  ri/R is the dimensionless radius-vector, 
which define i-th particle position; Γt  6 RTch/DV is 
the transition friction coefficient,  
 
 
2
5
1,
3 ( )N j j ij j ij ext
i
j j i ij
ρ u ρ u
h h  (4) 
 
is the magnetic field, which acts on the particle, 
ij  i – j; N is the number of the particles in the en-
semble; 3 4ext exth Η is the reduced external field; 
 T.V. LYUTYY, A.YU. POLYAKOV, S. DENISOV, P. HÄNGGI PROC. NAP 1, 04MFPN05 (2012) 
 
 
04MFPN05-2 
 
5
1,
7
( ) ( ) ( )
3
( )( )
15
N
ij j i i j ij j i ijdip
i
j j i ij
ij j ij i ij
ij
ρ u u u u ρ u u ρ
f
ρ u ρ u ρ
 (5) 
 
is the dipole force; and finally  
 
 
12 6
2
1,
24
N
ijsr
i
ij ijj j i ij
ρ
f  (6) 
 
is the force, which is modeling the repulsion. Here  
and  are phenomenological parameters, that define 
the depth of the potential well and distance between 
the particles centers, which corresponds to the minimal 
value of the Lennard-Jones potential, respectively. 
The random forces r ( ), and d, represent the parti-
cle interaction with heat bath, and are characterized by 
zero means and delta-like correlation functions 
 
 B( ) ( ) 3 ( ) ( , , )r r r chk T T , 
 B( ) ( ) 2 ( ) ( , , , )d d d chk T T x y z , 
 
where kB is the Boltzmann constant, T is the absolute 
temperature,  is Kronecker delta symbol, and ( ) is 
Dirac delta function. 
 
3. PARALLEL COMPUTATIONS WITH CUDA 
 
The most straightforward way to simulate ferrofluid 
dynamics is to use the so-called All-Pairs Algorithm, 
that takes into account interactions between all the 
pairs of the particles. The computation time in this 
case scales like O(N2),  thus making this brute-force 
approach not suitable for the numerical propagation of 
large systems, N >> 1.  
Computation time can be tangible decreased by im-
plementing the Barnes-Hut algorithm [4], which leads 
to the scaling of the computational time O(NlogN). The 
algorithm is based on the replacing of a group of parti-
cles by a pseudo-particle with integral characteristics 
of the group. If the group if far enough from the refer-
ence particle, then, instead of computing interaction 
between chosen particle and each particle in the group, 
it is reasonable to calculate the influence of the pseudo-
particle. The partition of the set of particles into the 
groups is performed in a hierarchical manner. The 
region where ferrofluid is located can be considered as 
a root cell, or a root group [4]. In 3d case the root cell is 
divided into 8 sub-cells, and each sub-cell is again di-
vided, etc., until every sub-sub-... -cell holds one or 
none particles. From every sub-cell we produce a pseu-
do particle has a magnetic moment u  equal to the sum 
of ui, that belong to the subcell. The position  of each 
pseudo particle determines as   i/N  (i  1,…,N ) 
Here N  is the number of particles in the subcell. 
The Barnes-Hut algorithm by the construction pos-
sesses a high-potential for parallelization [5] , but its 
implementation is quite challenging due to hierarchical 
tree it is based on. So we used recommendations in [5] 
to overcome these difficulties. 
Fig. 1 shows the efficiency of the Barnes-Hut GPU-
based algorithm. For the N  106 we observed a speed-
up of the order 6 105 as compared to the performance 
of the standard all-pairs algorithm on a CPU, and 
speed-up of the order 200 when compared to the all-
pairs implemented on a GPU. However, when the 
number of particles is relatively small, typically less 
than several thousands, the Barnes-Hut algorithm 
performs slower than the all-pairs  algorithm. This is 
because the tree-building procedure takes most of the 
computational time. 
 
 
 
 
Fig. 1 – Performance of different algorithms. 
 
4. TESTING OF THE DEVELOPED METHOD  
 
We performed simulations for the ensemble of ma-
ghemite ( -Fe2O3) particles,  by using Eqs.(1)-(3), with 
the parameters   3.1 105 A/m, D  5000 kg/m3, sus-
pended in a water of viscosity   0,89 10 – 3 Pa s.  
According to Shliomis [6], if the diluted (volume 
fraction  < 1 %) ferrofluid contains small particles for 
which the magnetic energy does not exceed the thermal 
energy, its magnetization is given by the Langevin 
function  
 
 0
B
coth 1 ,z
ext
zmH
k T
u . (7) 
 
The result of our simulation (see Fig. 2) are in a 
good agreement with Eq.(7), that is confirm the pro-
posed approach validity. The slight excess of the fer-
rofluid reduced magnetization over the Langevin func-
tion with  grows is explained by increasing of the 
dipole field component, that directed along the external 
field extzH . 
When the particle magnetic energy is larger then 
the thermal one, the ferrofluid response to the external 
field becomes more pronounced (see Fig. 3). At the 
same time, the role of interparticle interaction becomes 
more important. In particular, due to the antiferro-
magnetic coupling between the column clusters, and 
the formation of the ring-like clusters, the magnetiza-
tion of a dense ferrofluid can be several times smaller 
than the applied magnetic field. 
 
1.00E-01
1.00E+01
1.00E+03
1.00E+05
1.00E+07
1.00E+03 1.00E+04 1.00E+05 1.00E+06
All Pairs CPU
All Pairs GPU
Barnes-Hut GPU
N 
Time, ms 
 TECHNIQUE OF THE FAST FERROFLUIDS SIMULATION PROC. NAP 1, 04MFPN05 (2012) 
 
 
04MFPN05-3 
 
 
Fig. 2 – The magnetization curve for diluted (  < 1 %) ferro-
fluid. N  106, T  298 K  
 
5. CONCLUSIONS  
 
We proposed the new, GPU-based approach to fer-
rofluid simulations. By using the the Barnes-Hut algo-
rithm, which takes into account both the nearest-
neighbor correlation and interaction of the remote par-
ticles, we have reached an impressive calculation 
speed-up. One of the intriguing perspectives for further 
studies is to explore the dynamics of ferroliquids in 
thin vessels, a problem which is very important in the 
context of biomedical applications. 
 
 
Fig. 3 – The magnetization curves for different volume frac-
tion, when a > 1. N  106, T  298 K 
 
AKNOWLEDGEMENTS 
 
T.V.L. and A.Yu.P. acknowledge the support of the 
Cabinet of Ministers of Ukraine, The Program of 
Studying and Training for Students, PhD Students 
and Professor's Stuff Abroad (order dated April 13, 
2011 No 411) and the support of the Ministry of Edu-
cation, Science, Youth, and Sport of Ukraine (Project 
No 0112U001383).  
 
REFERENCES 
 
1. R. Rosensweig, Ferrohydrodynamics, (Cambridge Univer-
sity Press: 1985). 
2. Q. Pankhurst, J. Connolly, S. Jones, J. Dobson, J. Phys. 
D.: Appl. Phys. 36, R167 (2003). 
3. Z. Wang, C. Holm, H.W. Muller, Phys. Rev. E 66, 021405 
(2002). 
4. Josh Barnes, and Piet Hut, Nature 324, 446 (1986). 
5. Wen-mei W. Hwu, GPU Computing Gems Emerald Edi-
tion, 1st edition (Morgan Kaufmann Publishers Inc.: San 
Francisco, CA, USA: 2011). 
6. M. Shliomis, Sov. Phys. Usp. 17, 153 (1974). 
 
0
0.2
0.4
0.6
0.8
1
0 2 4 6 8
Simulation
Langevin function
0.00
0.20
0.40
0.60
0.80
1.00
0 0.05 0.1 0.15 0.2 0.25 0.3
Volume fraction 1%
Volume fraction 5%
Langevin function
uz  
 
 
uz  
 
uz  


Technique of the Fast Ferrofluids Simulation

Electronic Sumy State University Institutional Repository

PROCEEDINGS OF THE INTERNATIONAL CONFERENCE NANOMATERIALS: APPLICATIONS AND PROPERTIES Vol. 1 No 1, 04MFPN05(3pp) (2012)   2304-1862/2012/1(1)04MFPN05(3) 04MFPN05-1  2012 Sumy State University Technique of the Fast Ferrofluids Simulation  A.Yu. Polyakov1,2,*, T.V. Lyutyy1,†, S. Denisov2, P. Hänggi2  1 Sumy State University, 2, Rymsky Korsakov Str., 40007 Sumy, Ukraine 2 Institut für Physik, Universität Augsburg, Universitätsstraβe 1, D-86135 Augsburg, Germany  (Received 10 July 2012; published online 22 August 2012)  We present a highly-parallel implementation of the Langevin simulation method for modeling ferroflu-ids on Graphical Processor Units (GPU). Our method is based on the Barnes-Hut algorithm. As a bench-mark we use the straightforward 'all-to-all interaction' algorithm. The obtained results are in good agree-ment with known theoretical model. With the proposed method we were able to follow the evolution of a system of one million interacting particles over long time-scales, the task hitherto is out of reach with the standard, CPU-based numerical schemes.  Keywords: Ferrofluids, Molecular dynamics, N-body problem, Graphics processing unit, CUDA, Barnes-Hut algorithm.   PACS numbers: 47.65.Cb, 02.70.Ns                                                             * alexander.p1987@gmail.com † lyutyy@oeph.sumdu.edu.ua 1. INTRODUCTION  A system of N interacting particles is one of the most common set-ups in physics. Since analytical solu-tions typically do not exist when the number of inter-acted particles is larger then two, N > 2, the various numerical methods are used in order to get an insight in the dynamics of the corresponding system. The in-teraction involves all the pairs of particles, and there-fore the computational time and efforts scale like N2. This effect certainly constitutes a problem, and, in order to run big systems, researchers relied on big computers. Practically that means the use of large computational clusters consisting of many Central Processing Units (CPUs). Such clusters are physically big, expensive and consume a lot of electricity. Over the last decade, an alternative, Graphic Processing Units (GPUs), become more and more popular. GPUs, initial-ly designed to serve as the data-pipelines for graphical information, posses a high level of parallelism, while being cheap (as compared to CPUs) and small. The high computation capability of GPUs is owing to their architecture, which represents a set of hundreds of computational cores (note that the best modern CPUs have typically up to 16 cores), which can perform com-putations in parallel. Nowadays, the scientific GPU-computing is used in many areas of computational physics, thanks to the Compute United Device Archi-tecture (CUDA) developed by NVIDIA Corporation, which step simplified the development of programs for GPUs.  2. MODEL OF THE DISPERSED FERROFLUID  Theory of the continuous bulk ferrofluid is well-developed within the concept of ferrohydrodynamics [1]. However, recently the aspects of the dynamics of confined ferrofluids attracted a lot of attention because of the promising potential of ferrofluids for biotechnologies and medicine [2]. In this context, the disperse structure of a ferrofluidic media certainly deserves an attention. The dynamics of a ferromagnetic particle, suspended in liq-uid, can be described by a Langevin equation [3] that accounts for the long-range dipole interactions with other particles, caused by the particle magnetization, and for short-range steric repulsions, caused by the anti-agglomeration coating. Here we use the simplest model and consider an ensemble of identical spherical particles of radius R, submerged into the liquid of viscosity . We suppose also that the hydrodynamic radius is equal to R, and every particle has a homogenous structure with respect to its inertia and magnetic properties. Therefore, the dynamics of i-th particle is determined by the follow-ing equations   222 22,5i i chix iy iy ix r rTd du h u hd d VR D (1)  22222( )sin5( )cos ,iiy iz iz iy ki chiz ix ix iz k r rdu h u hdTdu h u hd VR D (2)   222 2,i idip chsr t di iTd dd d VR Dρ ρf f ξ  (3)  where i and  i are the azimuthal and polar angles of i-th particle;   t/Tch is the dimensionless time (Tch   R/ (3D/4 0)0.5;  is the particle material magnetization, 0  4 10 – 7H/m is the magnetic constant, D is the material particle density); ui  mi/m is the dimension-less magnetic moment (m  4 R3 /3  V , V is the par-ticle volume); Γr  8 RTch/DV is the rotation friction coefficient; i  ri/R is the dimensionless radius-vector, which define i-th particle position; Γt  6 RTch/DV is the transition friction coefficient,    251,3 ( )N j j ij j ij extij j i ijρ u ρ uh h  (4)  is the magnetic field, which acts on the particle, ij  i – j; N is the number of the particles in the en-semble; 3 4ext exth Η is the reduced external field;  T.V. LYUTYY, A.YU. POLYAKOV, S. DENISOV, P. HÄNGGI PROC. NAP 1, 04MFPN05 (2012)   04MFPN05-2  51,7( ) ( ) ( )3( )( )15Nij j i i j ij j i ijdipij j i ijij j ij i ijijρ u u u u ρ u u ρfρ u ρ u ρ (5)  is the dipole force; and finally    12 621,24Nijsriij ijj j i ijρf  (6)  is the force, which is modeling the repulsion. Here  and  are phenomenological parameters, that define the depth of the potential well and distance between the particles centers, which corresponds to the minimal value of the Lennard-Jones potential, respectively. The random forces r ( ), and d, represent the parti-cle interaction with heat bath, and are characterized by zero means and delta-like correlation functions   B( ) ( ) 3 ( ) ( , , )r r r chk T T ,  B( ) ( ) 2 ( ) ( , , , )d d d chk T T x y z ,  where kB is the Boltzmann constant, T is the absolute temperature,  is Kronecker delta symbol, and ( ) is Dirac delta function.  3. PARALLEL COMPUTATIONS WITH CUDA  The most straightforward way to simulate ferrofluid dynamics is to use the so-called All-Pairs Algorithm, that takes into account interactions between all the pairs of the particles. The computation time in this case scales like O(N2),  thus making this brute-force approach not suitable for the numerical propagation of large systems, N >> 1.  Computation time can be tangible decreased by im-plementing the Barnes-Hut algorithm [4], which leads to the scaling of the computational time O(NlogN). The algorithm is based on the replacing of a group of parti-cles by a pseudo-particle with integral characteristics of the group. If the group if far enough from the refer-ence particle, then, instead of computing interaction between chosen particle and each particle in the group, it is reasonable to calculate the influence of the pseudo-particle. The partition of the set of particles into the groups is performed in a hierarchical manner. The region where ferrofluid is located can be considered as a root cell, or a root group [4]. In 3d case the root cell is divided into 8 sub-cells, and each sub-cell is again di-vided, etc., until every sub-sub-... -cell holds one or none particles. From every sub-cell we produce a pseu-do particle has a magnetic moment u  equal to the sum of ui, that belong to the subcell. The position  of each pseudo particle determines as   i/N  (i  1,…,N ) Here N  is the number of particles in the subcell. The Barnes-Hut algorithm by the construction pos-sesses a high-potential for parallelization [5] , but its implementation is quite challenging due to hierarchical tree it is based on. So we used recommendations in [5] to overcome these difficulties. Fig. 1 shows the efficiency of the Barnes-Hut GPU-based algorithm. For the N  106 we observed a speed-up of the order 6 105 as compared to the performance of the standard all-pairs algorithm on a CPU, and speed-up of the order 200 when compared to the all-pairs implemented on a GPU. However, when the number of particles is relatively small, typically less than several thousands, the Barnes-Hut algorithm performs slower than the all-pairs  algorithm. This is because the tree-building procedure takes most of the computational time.     Fig. 1 – Performance of different algorithms.  4. TESTING OF THE DEVELOPED METHOD   We performed simulations for the ensemble of ma-ghemite ( -Fe2O3) particles,  by using Eqs.(1)-(3), with the parameters   3.1 105 A/m, D  5000 kg/m3, sus-pended in a water of viscosity   0,89 10 – 3 Pa s.  According to Shliomis [6], if the diluted (volume fraction  < 1 %) ferrofluid contains small particles for which the magnetic energy does not exceed the thermal energy, its magnetization is given by the Langevin function    0Bcoth 1 ,zextzmHk Tu . (7)  The result of our simulation (see Fig. 2) are in a good agreement with Eq.(7), that is confirm the pro-posed approach validity. The slight excess of the fer-rofluid reduced magnetization over the Langevin func-tion with  grows is explained by increasing of the dipole field component, that directed along the external field extzH . When the particle magnetic energy is larger then the thermal one, the ferrofluid response to the external field becomes more pronounced (see Fig. 3). At the same time, the role of interparticle interaction becomes more important. In particular, due to the antiferro-magnetic coupling between the column clusters, and the formation of the ring-like clusters, the magnetiza-tion of a dense ferrofluid can be several times smaller than the applied magnetic field.  1.00E-011.00E+011.00E+031.00E+051.00E+071.00E+03 1.00E+04 1.00E+05 1.00E+06All Pairs CPUAll Pairs GPUBarnes-Hut GPUN Time, ms  TECHNIQUE OF THE FAST FERROFLUIDS SIMULATION PROC. NAP 1, 04MFPN05 (2012)   04MFPN05-3   Fig. 2 – The magnetization curve for diluted (  < 1 %) ferro-fluid. N  106, T  298 K   5. CONCLUSIONS   We proposed the new, GPU-based approach to fer-rofluid simulations. By using the the Barnes-Hut algo-rithm, which takes into account both the nearest-neighbor correlation and interaction of the remote par-ticles, we have reached an impressive calculation speed-up. One of the intriguing perspectives for further studies is to explore the dynamics of ferroliquids in thin vessels, a problem which is very important in the context of biomedical applications.   Fig. 3 – The magnetization curves for different volume frac-tion, when a > 1. N  106, T  298 K  AKNOWLEDGEMENTS  T.V.L. and A.Yu.P. acknowledge the support of the Cabinet of Ministers of Ukraine, The Program of Studying and Training for Students, PhD Students and Professor's Stuff Abroad (order dated April 13, 2011 No 411) and the support of the Ministry of Edu-cation, Science, Youth, and Sport of Ukraine (Project No 0112U001383).   REFERENCES  1. R. Rosensweig, Ferrohydrodynamics, (Cambridge Univer-sity Press: 1985). 2. Q. Pankhurst, J. Connolly, S. Jones, J. Dobson, J. Phys. D.: Appl. Phys. 36, R167 (2003). 3. Z. Wang, C. Holm, H.W. Muller, Phys. Rev. E 66, 021405 (2002). 4. Josh Barnes, and Piet Hut, Nature 324, 446 (1986). 5. Wen-mei W. Hwu, GPU Computing Gems Emerald Edi-tion, 1st edition (Morgan Kaufmann Publishers Inc.: San Francisco, CA, USA: 2011). 6. M. Shliomis, Sov. Phys. Usp. 17, 153 (1974).  00.20.40.60.810 2 4 6 8SimulationLangevin function0.000.200.400.600.801.000 0.05 0.1 0.15 0.2 0.25 0.3Volume fraction 1%Volume fraction 5%Langevin functionuz    uz   uz  

https://essuir.sumdu.edu.ua/bitstream/123456789/35341/1/princon_2012_1_4_22.pdf

Technique of the Fast Ferrofluids Simulation

Abstract

Similar works

Full text

Available Versions

SocioEconomic Challenges Journal (Sumy State University)

Electronic Sumy State University Institutional Repository