Commit 6ab6af86 authored by Dr. Bogdan Tanygin's avatar Dr. Bogdan Tanygin Committed by GitHub

Merge pull request #141 from espressomd/python

Lb force fix (#1099)
parents ce730fe6 c981453d
......@@ -76,7 +76,7 @@ for i in range(40):
print("\nIntegration finished.")
# get force that is exerted on the sphere
force = sphere.get_params()["force"]
force = sphere.get_force()
print("Measured force: f=%f" %size(force))
stokes_force = 6*np.pi*kinematic_visc*radius*size(v)
......
......@@ -3,6 +3,7 @@
#include <memory>
#include "config.hpp"
#include "shapes/NoWhere.hpp"
#include "shapes/Shape.hpp"
#include "Vector.hpp"
......
......@@ -25,3 +25,4 @@ if any(i in espressomd.code_info.features() for i in ["LB_BOUNDARIES", "LB_BOUND
class LBBoundary(ScriptInterfaceHelper):
_so_name = "LBBoundaries::LBBoundary"
_so_bind_methods = ("get_force",)
......@@ -12,7 +12,6 @@ class LBBoundary : public AutoParameters {
public:
LBBoundary() : m_lbboundary(new ::LBBoundaries::LBBoundary()) {
add_parameters({{"velocity", m_lbboundary->velocity()},
{"force", m_lbboundary->get_force()},
{"shape",
[this](Variant const &value) {
m_shape =
......@@ -26,21 +25,28 @@ public:
return (m_shape != nullptr) ? m_shape->id() : ObjectId();
}}});
#ifdef EK_BOUNDARIES
add_parameters({
{"charge_density",
[this](Variant const &value) {
m_lbboundary->set_charge_density(boost::get<double>(value));
},
[this]() { return m_lbboundary->charge_density(); }},
{"net_charge",
[this](Variant const &value) {
m_lbboundary->set_net_charge(boost::get<double>(value));
},
[this]() { return m_lbboundary->net_charge(); }},
});
add_parameters(
{{"charge_density",
[this](Variant const &value) {
m_lbboundary->set_charge_density(boost::get<double>(value));
},
[this]() { return m_lbboundary->charge_density(); },
VariantType::DOUBLE, 0},
{"net_charge",
[this](Variant const &value) {
m_lbboundary->set_net_charge(boost::get<double>(value));
},
[this]() { return m_lbboundary->net_charge(); }, VariantType::DOUBLE,
0}});
#endif
}
Variant call_method(const std::string &name, const VariantMap &) {
if (name == "get_force") {
return m_lbboundary->get_force();
}
}
const std::string name() const override { return "LBBoundaries:LBBoundary"; }
std::shared_ptr<::LBBoundaries::LBBoundary> lbboundary() {
......
......@@ -17,6 +17,7 @@ set(py_tests bondedInteractions.py
rotational_inertia.py
lbgpu_remove_total_momentum.py
tabulated.py
lb_stokes_sphere_gpu.py
)
if(PY_H5PY)
set(py_tests ${py_tests} h5md.py)
......
### Measuring the force on a single sphere immersed in a fluid with
### fixed velocity boundary conditions created by two
### walls at finite distance.
### The force is compared to th analytical result F=6 pi eta r v
### i.e. the stokes force on the particles.
# We create a box of size box_width x box_width x box_length and
# place an object in the center. We measure the drag force
# in z direction. We create walls in the xz and yz plane at the box
# boundaries, where the velocity is fixed to $v.
#
from espressomd import System, lb, lbboundaries, shapes, has_features
import unittest as ut
import numpy as np
import sys
@ut.skipIf(not has_features(["LB_GPU", "LB_BOUNDARIES_GPU"]),
"Features not available, skipping test!")
class Stokes(ut.TestCase):
es = System()
def test_stokes(self):
# System setup
#system = System()
system=self.es
agrid = 1
radius = 5.5
box_width = 64
real_width = box_width+2*agrid
box_length = 64
system.box_l = [real_width, real_width, box_length]
system.time_step = 0.2
system.cell_system.skin = 0.4
# The temperature is zero.
system.thermostat.set_lb(kT=0)
# LB Parameters
v = [0,0,0.01] # The boundary slip
kinematic_visc = 1.0
# Invoke LB fluid
lbf = lb.LBFluid_GPU(visc=kinematic_visc, dens=1, agrid=agrid, tau=system.time_step, fric=1)
system.actors.add(lbf)
# Setup walls
walls = [None] * 4
walls[0] = lbboundaries.LBBoundary(shape=shapes.Wall(normal=[-1,0,0],
dist = -(1+box_width)), velocity=v)
walls[1] = lbboundaries.LBBoundary(shape=shapes.Wall(normal=[1,0,0], dist = 1),
velocity=v)
walls[2] = lbboundaries.LBBoundary(shape=shapes.Wall(normal=[0,-1,0],
dist = -(1+box_width)), velocity=v)
walls[3] = lbboundaries.LBBoundary(shape=shapes.Wall(normal=[0,1,0], dist = 1),
velocity=v)
for wall in walls:
system.lbboundaries.add(wall)
# setup sphere without slip in the middle
sphere = lbboundaries.LBBoundary(shape=shapes.Sphere(radius=radius,
center = [real_width/2] * 2 + [box_length/2],
direction = 1))
system.lbboundaries.add(sphere)
def size(vector):
tmp = 0
for k in vector:
tmp+=k*k
return np.sqrt(tmp)
system.integrator.run(4000)
print("\nIntegration finished.")
# get force that is exerted on the sphere
force = sphere.get_force()
print("Measured force: f=%f" %size(force))
stokes_force = 6*np.pi*kinematic_visc*radius*size(v)
print("Stokes' Law says: f=%f" %stokes_force)
self.assertTrue(abs(1.0 - size(force)/stokes_force) <0.06)
if __name__ == "__main__":
ut.main()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment