RYUON-stokes
Stokesian dynamics simulator

RYUON| documents| stokes| libstokes libiter| twobody
outline| visualization| download| wrappers| 1st

Outline

RYUON-stokes is a Stokesian dynamics simulator. Using libstokes, RYUON-stokes can solve resistance and mobility problems under the periodic boundary condition with F, FT, and FTS versions for full 3D as well as semi-2D (monolayer) configurations. For each of them, there are two functions, with or without lubrication correction. For mobility problem, also fixed particles can be handled.

Visualization

by POVRAY VTK Postscript

Dancing Particles in Shear Flow

6 particles on circle 10 particles in line

visualized by stnc2pov.py and POVRAY.

Download

please visit File Releases @ SF.

Changes and Plans

Requirements

0) using libstokes from scripting languages

You can use libstokes from many programming languages. (All codes are in stokes-0.4.tar.bz2.)

1) the first problem

Let's solve a resistance problem of 8 particles with the same velocity placed at SC lattice cites. Of course, the expected forces would be same for all particles. (All codes are in stokes-0.4.tar.bz2.)

Description C (html) Python (html) Ruby (html) Perl (html) Guile (html) Octave

preamble

 
#include <stdio.h>
#include <stdlib.h>

#include <libstokes.h>
#include <libiter.h>

/* main program */
int
main (int argc, char** argv)
{
 
import stokes
 
require 'stokes'
 
use stokes;
 
(load-extension "./stokes.so" "SWIG_init")

initialize struct stokes * sys.

 
  struct stokes *sys = stokes_init ();
 
sys = stokes.stokes_init()
 
sys = Stokes::stokes_init()
 
$sys = stokes::stokes_init();
 
(define sys (stokes-init))

set np and nm


(you must call stokes_set_np() because this also allocate the memory for pos.) 
  int np, nm;
  np = 8;
  nm = 8;
  stokes_set_np (sys, np, nm);
 
np = 8
nm = 8
stokes.stokes_set_np(sys, np, nm)
 
np = 8
nm = 8
Stokes::stokes_set_np(sys, np, nm)
# you must call stokes_set_np() because
# this also allocate the memory for pos.
 
$np = 8;
$nm = 8;
stokes::stokes_set_np($sys, $np, $nm);
# you must call stokes_set_np() because
# this also allocate the memory for pos.
 
(define np 8)
(define nm 8)
(stokes-set-np sys np nm)

set periodic lattice parameters.

 
  double lx, ly, lz;
  lx = 10.0;
  ly = 10.0;
  lz = 10.0;
  stokes_set_l (sys, lx, ly, lz);

  tratio = 60.25;
  xi = xi_by_tratio (sys, tratio);

  cutlim = 1.0e-12;
  stokes_set_xi (sys, xi, cutlim);

  fprintf (stdout, "xi = %f\n", xi);

  sys->lubcut = 2.0000000001;
 
lx = 10.0
ly = 10.0
lz = 10.0
stokes.stokes_set_l(sys, lx, ly, lz)

tratio = 60.25
xi = stokes.xi_by_tratio(sys, tratio)
cutlim = 1.0e-12
stokes.stokes_set_xi(sys, xi, cutlim)

print 'xi =', xi

sys.lubcut = 2.0000000001
 
lx = 10.0
ly = 10.0
lz = 10.0
Stokes::stokes_set_l(sys, lx, ly, lz)
# you must call stokes_set_l() becuase
# this also initialize parameters other than lx,ly,lz.

tratio = 60.25
xi = Stokes::xi_by_tratio(sys, tratio)

cutlim = 1.0e-12
Stokes::stokes_set_xi(sys, xi, cutlim)

print "xi = ", xi, "\n"

sys.lubcut = 2.0000000001
 
$lx = 10.0;
$ly = 10.0;
$lz = 10.0;
stokes::stokes_set_l($sys, $lx, $ly, $lz);

$tratio = 60.25;
$xi = stokes::xi_by_tratio($sys, $tratio);
$cutlim = 1.0e-12;
stokes::stokes_set_xi($sys, $xi, $cutlim);

print "xi = ", $xi, "\n";

$sys->{lubcut} = 2.0000000001;
 
(define lx 10.0)
(define ly 10.0)
(define lz 10.0)
(stokes-set-l sys lx ly lz)

(define tratio 60.25)
(define xi (xi-by-tratio sys tratio))

(define cutlim 1.0e-12)
(stokes-set-xi sys xi cutlim)

(display "xi = ")
(display xi)
(newline)

;(stokes-lubcut-set sys 2.0000000001)
;(stokes-lubcut-get sys)
; with '-emit-setter' you can write those as
(set! (stokes-lubcut sys) 2.0000000001)
(stokes-lubcut sys)
 
l = [10.0, 10.0, 10.0]

set the iterative solver.

 
  stokes_set_iter (sys, "gmres", 2000, 20, 1.0e-6, 1, stdout);
 
stokes.stokes_set_iter(sys, "gmres", 2000, 20, 1.0e-6,
                       1, stokes.get_stdout())
 
Stokes::stokes_set_iter(sys, "gmres", 2000, 20, 1.0e-6,
                        1, Stokes::get_stdout())
 
stokes::stokes_set_iter($sys, "gmres", 2000, 20, 1.0e-6,
			1, stokes::get_stdout());
 
(stokes-set-iter sys "gmres" 2000 20 1.0e-6 1 (get-stdout))

allocate memories for pos, u, and f and initialize them.

 
  int i;
  double * pos;
  double * u;
  double * f;
  pos = (double *) calloc (np * 3, sizeof (double));
  u   = (double *) calloc (np * 3, sizeof (double));
  f   = (double *) calloc (np * 3, sizeof (double));

  pos[ 0] = 0.0; // x component
  pos[ 1] = 0.0; // y component
  pos[ 2] = 0.0; // z component
  pos[ 3] = 5.0; pos[ 4] = 0.0; pos[ 5] = 0.0;
  pos[ 6] = 0.0; pos[ 7] = 5.0; pos[ 8] = 0.0;
  pos[ 9] = 0.0; pos[10] = 0.0; pos[11] = 5.0;
  pos[12] = 5.0; pos[13] = 5.0; pos[14] = 0.0;
  pos[15] = 0.0; pos[16] = 5.0; pos[17] = 5.0;
  pos[18] = 5.0; pos[19] = 0.0; pos[20] = 5.0;
  pos[21] = 5.0; pos[22] = 5.0; pos[23] = 5.0;

  for (i = 0; i < np * 3; i ++)
    {
      u[i] = 1.0;
    }
 
pos = stokes.darray(np*3)
u   = stokes.darray(np*3)
f   = stokes.darray(np*3)

pos[ 0] = 0.0 # x component
pos[ 1] = 0.0 # y component
pos[ 2] = 0.0 # z component

pos[ 3] = 5.0
pos[ 4] = 0.0
pos[ 5] = 0.0

pos[ 6] = 0.0
pos[ 7] = 5.0
pos[ 8] = 0.0

pos[ 9] = 0.0
pos[10] = 0.0
pos[11] = 5.0

pos[12] = 5.0
pos[13] = 5.0
pos[14] = 0.0

pos[15] = 0.0
pos[16] = 5.0
pos[17] = 5.0

pos[18] = 5.0
pos[19] = 0.0
pos[20] = 5.0

pos[21] = 5.0
pos[22] = 5.0
pos[23] = 5.0

for i in range(np*3):
    u[i] = 1.0
 
pos = Stokes::Darray.new(np*3)
u   = Stokes::Darray.new(np*3)
f   = Stokes::Darray.new(np*3)

# set pos in SC lattice configuration
pos[ 0] = 0.0 # x component
pos[ 1] = 0.0 # y component
pos[ 2] = 0.0 # z component

pos[ 3] = 5.0
pos[ 4] = 0.0
pos[ 5] = 0.0

pos[ 6] = 0.0
pos[ 7] = 5.0
pos[ 8] = 0.0

pos[ 9] = 0.0
pos[10] = 0.0
pos[11] = 5.0

pos[12] = 5.0
pos[13] = 5.0
pos[14] = 0.0

pos[15] = 0.0
pos[16] = 5.0
pos[17] = 5.0

pos[18] = 5.0
pos[19] = 0.0
pos[20] = 5.0

pos[21] = 5.0
pos[22] = 5.0
pos[23] = 5.0

for i in 0..(np*3-1)
    u[i] = 1.0
end
 
$pos = new stokes::darray($np*3);
$u   = new stokes::darray($np*3);
$f   = new stokes::darray($np*3);

$pos->setitem( 0, 0.0); # x component
$pos->setitem( 1, 0.0); # y component
$pos->setitem( 2, 0.0); # z component
$pos->setitem( 3, 5.0); $pos->setitem( 4, 0.0); $pos->setitem( 5, 0.0);
$pos->setitem( 6, 0.0); $pos->setitem( 7, 5.0); $pos->setitem( 8, 0.0);
$pos->setitem( 9, 0.0); $pos->setitem(10, 0.0); $pos->setitem(11, 5.0);
$pos->setitem(12, 5.0); $pos->setitem(13, 5.0); $pos->setitem(14, 0.0);
$pos->setitem(15, 0.0); $pos->setitem(16, 5.0); $pos->setitem(17, 5.0);
$pos->setitem(18, 5.0); $pos->setitem(19, 0.0); $pos->setitem(20, 5.0);
$pos->setitem(21, 5.0); $pos->setitem(22, 5.0); $pos->setitem(23, 5.0);

for ($i=0; $i < $np*3; $i++){
    $u->setitem($i, 1.0);
}
 
(define pos (new-darray (* np 3)))
(define u   (new-darray (* np 3)))
(define f   (new-darray (* np 3)))

; set pos in SC lattice configuration
(darray-setitem pos 0  0.0) ; x component
(darray-setitem pos 1  0.0) ; y component
(darray-setitem pos 2  0.0) ; z component

(darray-setitem pos 3  5.0)
(darray-setitem pos 4  0.0)
(darray-setitem pos 5  0.0)

(darray-setitem pos 6  0.0)
(darray-setitem pos 7  5.0)
(darray-setitem pos 8  0.0)

(darray-setitem pos 9  0.0)
(darray-setitem pos 10 0.0)
(darray-setitem pos 11 5.0)

(darray-setitem pos 12 5.0)
(darray-setitem pos 13 5.0)
(darray-setitem pos 14 0.0)

(darray-setitem pos 15 0.0)
(darray-setitem pos 16 5.0)
(darray-setitem pos 17 5.0)

(darray-setitem pos 18 5.0)
(darray-setitem pos 19 0.0)
(darray-setitem pos 20 5.0)

(darray-setitem pos 21 5.0)
(darray-setitem pos 22 5.0)
(darray-setitem pos 23 5.0)

; set u and f
(do ((i 0 (1+ i)))
    ((>= i (* np 3)))
  (darray-setitem u i 1.0))
 
pos = [0.0, 0.0, 0.0,\
	0.0, 0.0, 5.0,\
	0.0, 5.0, 0.0,\
	0.0, 5.0, 5.0,\
	5.0, 0.0, 0.0,\
	5.0, 0.0, 5.0,\
	5.0, 5.0, 0.0,\
	5.0, 5.0, 5.0]
u = [1.0, 1.0, 1.0,\
      1.0, 1.0, 1.0,\
      1.0, 1.0, 1.0,\
      1.0, 1.0, 1.0,\
      1.0, 1.0, 1.0,\
      1.0, 1.0, 1.0,\
      1.0, 1.0, 1.0,\
      1.0, 1.0, 1.0]

print pos, u, and f for check.

 
  fprintf (stdout, "pos:\n");
  for (i = 0; i < np; i ++)
    {
      fprintf (stdout, "%d %f %f %f\n",
	       i,
	       pos[i*3 + 0],
	       pos[i*3 + 1],
	       pos[i*3 + 2]);
    }

  fprintf (stdout, "u:\n");
  for (i = 0; i < np; i ++)
    {
      fprintf (stdout, "%d %f %f %f\n",
	       i,
	       u[i*3 + 0],
	       u[i*3 + 1],
	       u[i*3 + 2]);
    }
 
print 'pos:'
for i in range(np):
    print i,pos[i*3],pos[i*3+1],pos[i*3+2]

print 'u:'
for i in range(np):
    print i, u[i*3], u[i*3+1], u[i*3+2]
 
print "pos:\n"
for i in 0..(np-1)
    print i, ' ', pos[i*3], ' ', pos[i*3+1], ' ', pos[i*3+2], "\n"
end

print "u:\n"
for i in 0..(np-1)
    print i, ' ', u[i*3], ' ', u[i*3+1], ' ', u[i*3+2], "\n"
end
 
print "pos:\n";
for ($i=0; $i < $np; $i++){
    printf ("%d %f %f %f\n",
	    $i,
	    $pos->getitem($i*3),
	    $pos->getitem($i*3+1),
	    $pos->getitem($i*3+2));
}
print "u:\n";

for ($i=0; $i < $np; $i++){
    printf ("%d %f %f %f\n",
	    $i,
	    $u->getitem($i*3),
	    $u->getitem($i*3+1),
	    $u->getitem($i*3+2));
}
 
(display "pos:")
(newline)
(display-darray3 pos np)

(display "u:")
(newline)
(display-darray3 u np)

call the solver solve_res_ewald_3f(). (the real calculation is done here.)

 
  stokes_set_pos (sys, pos);
  solve_res_ewald_3f (sys, u, f);
 
stokes.stokes_set_pos(sys, pos)
stokes.solve_res_ewald_3f(sys, u, f)
 
Stokes::stokes_set_pos(sys, pos)
Stokes::calc_res_ewald_3f(sys, u, f)
 
stokes::stokes_set_pos($sys, $pos);
stokes::solve_res_ewald_3f($sys, $u, $f);
 
(set! (stokes-pos sys) pos)

(solve-res-ewald-3f sys u f)
 
f = stokes_res_ewald_3f(pos, u, l, 60.25)

print the result f.

 
  fprintf (stdout, "f:\n");
  for (i = 0; i < np; i ++)
    {
      fprintf (stdout, "%d %f %f %f\n",
	       i,
	       f[i*3 + 0],
	       f[i*3 + 1],
	       f[i*3 + 2]);
    }

  return 0;
}
 
print 'f:'
for i in range(np):
    print i, f[i*3], f[i*3+1], f[i*3+2]
 
print "f:\n"
for i in 0..(np-1)
    print i, ' ', f[i*3], ' ', f[i*3+1], ' ', f[i*3+2], "\n"
end
 
print "f:\n";
for ($i=0; $i < $np; $i++){
    printf ("%d %f %f %f\n",
	    $i,
	    $f->getitem($i*3),
	    $f->getitem($i*3+1),
	    $f->getitem($i*3+2));
}
 
(display "f:")
(newline)
(display-darray3 f np)

Results:

 
xi = 0.350941
pos:
0 0.000000 0.000000 0.000000
1 5.000000 0.000000 0.000000
2 0.000000 5.000000 0.000000
3 0.000000 0.000000 5.000000
4 5.000000 5.000000 0.000000
5 0.000000 5.000000 5.000000
6 5.000000 0.000000 5.000000
7 5.000000 5.000000 5.000000
u:
0 1.000000 1.000000 1.000000
1 1.000000 1.000000 1.000000
2 1.000000 1.000000 1.000000
3 1.000000 1.000000 1.000000
4 1.000000 1.000000 1.000000
5 1.000000 1.000000 1.000000
6 1.000000 1.000000 1.000000
7 1.000000 1.000000 1.000000
libiter: iter=20 res=2.799607e-16
f:
0 2.145689 2.145689 2.145689
1 2.145689 2.145689 2.145689
2 2.145689 2.145689 2.145689
3 2.145689 2.145689 2.145689
4 2.145689 2.145689 2.145689
5 2.145689 2.145689 2.145689
6 2.145689 2.145689 2.145689
7 2.145689 2.145689 2.145689
 
xi = 0.350941271209
pos:
0 0.0 0.0 0.0
1 5.0 0.0 0.0
2 0.0 5.0 0.0
3 0.0 0.0 5.0
4 5.0 5.0 0.0
5 0.0 5.0 5.0
6 5.0 0.0 5.0
7 5.0 5.0 5.0
u:
0 1.0 1.0 1.0
1 1.0 1.0 1.0
2 1.0 1.0 1.0
3 1.0 1.0 1.0
4 1.0 1.0 1.0
5 1.0 1.0 1.0
6 1.0 1.0 1.0
7 1.0 1.0 1.0
libiter: iter=20 res=2.799607e-16
f:
0 2.14568872011 2.14568872011 2.14568872011
1 2.14568872011 2.14568872011 2.14568872011
2 2.14568872011 2.14568872011 2.14568872011
3 2.14568872011 2.14568872011 2.14568872011
4 2.14568872011 2.14568872011 2.14568872011
5 2.14568872011 2.14568872011 2.14568872011
6 2.14568872011 2.14568872011 2.14568872011
7 2.14568872011 2.14568872011 2.14568872011
 
xi = 0.350941271209315
pos:
0 0.0 0.0 0.0
1 5.0 0.0 0.0
2 0.0 5.0 0.0
3 0.0 0.0 5.0
4 5.0 5.0 0.0
5 0.0 5.0 5.0
6 5.0 0.0 5.0
7 5.0 5.0 5.0
u:
0 1.0 1.0 1.0
1 1.0 1.0 1.0
2 1.0 1.0 1.0
3 1.0 1.0 1.0
4 1.0 1.0 1.0
5 1.0 1.0 1.0
6 1.0 1.0 1.0
7 1.0 1.0 1.0
libiter: iter=20 res=2.799607e-16
f:
0 2.14568872010948 2.14568872010949 2.14568872010948
1 2.14568872010948 2.14568872010949 2.14568872010949
2 2.14568872010948 2.14568872010949 2.14568872010948
3 2.14568872010948 2.14568872010949 2.14568872010949
4 2.14568872010948 2.14568872010949 2.14568872010949
5 2.14568872010948 2.14568872010949 2.14568872010948
6 2.14568872010948 2.14568872010949 2.14568872010949
7 2.14568872010948 2.14568872010949 2.14568872010949
 
xi = 0.350941271209315
pos:
0 0.000000 0.000000 0.000000
1 5.000000 0.000000 0.000000
2 0.000000 5.000000 0.000000
3 0.000000 0.000000 5.000000
4 5.000000 5.000000 0.000000
5 0.000000 5.000000 5.000000
6 5.000000 0.000000 5.000000
7 5.000000 5.000000 5.000000
u:
0 1.000000 1.000000 1.000000
1 1.000000 1.000000 1.000000
2 1.000000 1.000000 1.000000
3 1.000000 1.000000 1.000000
4 1.000000 1.000000 1.000000
5 1.000000 1.000000 1.000000
6 1.000000 1.000000 1.000000
7 1.000000 1.000000 1.000000
libiter: iter=20 res=2.799607e-16
f:
0 2.145689 2.145689 2.145689
1 2.145689 2.145689 2.145689
2 2.145689 2.145689 2.145689
3 2.145689 2.145689 2.145689
4 2.145689 2.145689 2.145689
5 2.145689 2.145689 2.145689
6 2.145689 2.145689 2.145689
7 2.145689 2.145689 2.145689
 
xi = 0.350941271209315
pos:
0 0.0 0.0 0.0
1 5.0 0.0 0.0
2 0.0 5.0 0.0
3 0.0 0.0 5.0
4 5.0 5.0 0.0
5 0.0 5.0 5.0
6 5.0 0.0 5.0
7 5.0 5.0 5.0
u:
0 1.0 1.0 1.0
1 1.0 1.0 1.0
2 1.0 1.0 1.0
3 1.0 1.0 1.0
4 1.0 1.0 1.0
5 1.0 1.0 1.0
6 1.0 1.0 1.0
7 1.0 1.0 1.0
libiter: iter=20 res=2.799607e-16
f:
0 2.14568872010948 2.14568872010949 2.14568872010948
1 2.14568872010948 2.14568872010949 2.14568872010949
2 2.14568872010948 2.14568872010949 2.14568872010948
3 2.14568872010948 2.14568872010949 2.14568872010949
4 2.14568872010948 2.14568872010949 2.14568872010949
5 2.14568872010948 2.14568872010949 2.14568872010948
6 2.14568872010948 2.14568872010949 2.14568872010949
7 2.14568872010948 2.14568872010949 2.14568872010949
 
pos =

  0  0  0  0  0  5  0  5  0  0  5  5  5  0  0  5  0  5  5  5  0  5  5  5

u =

  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1

l =

  10  10  10

libiter: iter=20 res=3.325435e-16
f =

  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457
  2.1457

syntax highlighted by Code2HTML, v. 0.9.1.

outline| visualization| download| wrappers| 1st
RYUON| documents| stokes| libstokes libiter| twobody

e-mail: kengoichiki@gmail.com
Copyright (C) 2006-2017 Kengo Ichiki. All rights reserved.