diff --git a/src/config/config.hpp b/src/config/config.hpp index c14042753b2..65ba5cae25a 100644 --- a/src/config/config.hpp +++ b/src/config/config.hpp @@ -100,9 +100,4 @@ #define MAX_OBJECTS_IN_FLUID 10000 #endif -/** Maximal number of soft particles per immersed boundary */ -#ifndef IBM_MAX_NUM -#define IBM_MAX_NUM 1000 -#endif - #endif diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index c31c472edac..cd1026ac61d 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -257,6 +257,12 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, if (auto const *iap = boost::get(&iaparams)) { return iap->energy(p1.r.p, p2->r.p, p3->r.p); } + if (auto const *iap = boost::get(&iaparams)) { + runtimeWarningMsg() << "Unsupported bond type " + + std::to_string(iaparams.which()) + + " in energy calculation."; + return 0.; + } throw BondUnknownTypeError(); } // 2 partners if (n_partners == 3) { @@ -266,6 +272,12 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, if (auto const *iap = boost::get(&iaparams)) { return iap->energy(p2->r.p, p1.r.p, p3->r.p, p4->r.p); } + if (auto const *iap = boost::get(&iaparams)) { + runtimeWarningMsg() << "Unsupported bond type " + + std::to_string(iaparams.which()) + + " in energy calculation."; + return 0.; + } throw BondUnknownTypeError(); } // 3 partners if (n_partners == 0) { diff --git a/src/core/immersed_boundary/ImmersedBoundaries.cpp b/src/core/immersed_boundary/ImmersedBoundaries.cpp index 081aa7e39d8..164e4392f95 100644 --- a/src/core/immersed_boundary/ImmersedBoundaries.cpp +++ b/src/core/immersed_boundary/ImmersedBoundaries.cpp @@ -32,10 +32,11 @@ #include #include -#include +#include -#include +#include #include +#include /** Calculate volumes, volume force and add it to each virtual particle. */ void ImmersedBoundaries::volume_conservation(CellStructure &cs) { @@ -69,7 +70,7 @@ void ImmersedBoundaries::init_volume_conservation(CellStructure &cs) { // accidentally during the integration. Then we must not reset the // reference BoundariesFound = true; - if (v->volRef == 0) { + if (v->volRef == 0.) { v->volRef = VolumesCurrent[v->softID]; } } @@ -98,7 +99,7 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) { return; // Partial volumes for each soft particle, to be summed up - std::vector tempVol(IBM_MAX_NUM); + std::vector tempVol(VolumesCurrent.size()); // Loop over all particles on local node cs.bond_loop( @@ -143,12 +144,10 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) { return false; }); - for (int i = 0; i < IBM_MAX_NUM; i++) - VolumesCurrent[i] = 0; - // Sum up and communicate - MPI_Allreduce(tempVol.data(), VolumesCurrent.data(), IBM_MAX_NUM, MPI_DOUBLE, - MPI_SUM, comm_cart); + boost::mpi::all_reduce(comm_cart, tempVol.data(), + static_cast(tempVol.size()), + VolumesCurrent.data(), std::plus()); } /** Calculate and add the volume force to each node */ diff --git a/src/core/immersed_boundary/ImmersedBoundaries.hpp b/src/core/immersed_boundary/ImmersedBoundaries.hpp index ba6d815b0a4..05f9250bdd0 100644 --- a/src/core/immersed_boundary/ImmersedBoundaries.hpp +++ b/src/core/immersed_boundary/ImmersedBoundaries.hpp @@ -23,15 +23,29 @@ #include "CellStructure.hpp" +#include +#include #include class ImmersedBoundaries { public: ImmersedBoundaries() : VolumeInitDone(false), BoundariesFound(false) { - VolumesCurrent.resize(IBM_MAX_NUM); + VolumesCurrent.resize(10); } void init_volume_conservation(CellStructure &cs); void volume_conservation(CellStructure &cs); + void register_softID(int softID) { + assert(softID >= 0); + auto const new_size = static_cast(softID) + 1; + if (new_size > VolumesCurrent.size()) { + VolumesCurrent.resize(new_size); + } + } + double get_current_volume(int softID) const { + assert(softID >= 0); + assert(softID < VolumesCurrent.size()); + return VolumesCurrent[softID]; + } private: void calc_volumes(CellStructure &cs); diff --git a/src/core/immersed_boundary/ibm_volcons.cpp b/src/core/immersed_boundary/ibm_volcons.cpp index 4ff52ad70c3..37415fbb55d 100644 --- a/src/core/immersed_boundary/ibm_volcons.cpp +++ b/src/core/immersed_boundary/ibm_volcons.cpp @@ -18,25 +18,31 @@ */ #include "ibm_volcons.hpp" -#include "errorhandling.hpp" + +#include "communication.hpp" +#include "immersed_boundaries.hpp" + +#include +#include + +void mpi_set_n_ibm_volcons_bonds_local(int softID) { + immersed_boundaries.register_softID(softID); +} + +REGISTER_CALLBACK(mpi_set_n_ibm_volcons_bonds_local) + +void mpi_set_n_ibm_volcons_bonds(int softID) { + mpi_call_all(mpi_set_n_ibm_volcons_bonds_local, softID); +} /** Set parameters of volume conservation */ IBMVolCons::IBMVolCons(const int softID, const double kappaV) { - // Specific stuff - if (softID > IBM_MAX_NUM) { - runtimeErrorMsg() << "Error: softID (" << softID - << ") is larger than IBM_MAX_NUM (" << IBM_MAX_NUM << ")"; - } - if (softID < 0) { - runtimeErrorMsg() << "Error: softID (" << softID - << ") must be non-negative"; - } - this->softID = softID; this->kappaV = kappaV; - volRef = 0; // NOTE: We cannot compute the reference volume here because not all // interactions are setup and thus we do not know which triangles belong to - // this softID. Calculate it later in the init function of \ref - // ImmersedBoundaries + // this softID. Calculate it later in the init function of + // \ref ImmersedBoundaries::init_volume_conservation() + volRef = 0.; + mpi_set_n_ibm_volcons_bonds(softID); } diff --git a/src/python/espressomd/interactions.pxd b/src/python/espressomd/interactions.pxd index 7d042dd950d..5ed7b3f42e4 100644 --- a/src/python/espressomd/interactions.pxd +++ b/src/python/espressomd/interactions.pxd @@ -528,7 +528,7 @@ cdef extern from "thermostat.hpp": cdef extern from "immersed_boundary/ImmersedBoundaries.hpp": cppclass ImmersedBoundaries: - pass + double get_current_volume(int softID) cdef extern from "immersed_boundaries.hpp": extern ImmersedBoundaries immersed_boundaries diff --git a/src/python/espressomd/interactions.pyx b/src/python/espressomd/interactions.pyx index ee082e6c19a..5acd4d3a86c 100644 --- a/src/python/espressomd/interactions.pyx +++ b/src/python/espressomd/interactions.pyx @@ -21,7 +21,7 @@ from cython.operator cimport dereference import collections include "myconfig.pxi" -from .utils import requires_experimental_features, is_valid_type +from .utils import is_valid_type from .utils cimport check_type_or_throw_except @@ -2942,7 +2942,6 @@ class IBM_Tribend(BondedInteraction): self._params["ind4"], self._params["kb"], flat)) -@requires_experimental_features("No test coverage") class IBM_VolCons(BondedInteraction): """ @@ -2954,7 +2953,8 @@ class IBM_VolCons(BondedInteraction): ---------- softID : :obj:`int` Used to identify the object to which this bond belongs. Each object - (cell) needs its own ID + (cell) needs its own ID. For performance reasons, it is best to + start from ``softID=0`` and increment by 1 for each subsequent bond. kappaV : :obj:`float` Modulus for volume force @@ -2988,6 +2988,14 @@ class IBM_VolCons(BondedInteraction): set_bonded_ia_params(self._bond_id, CoreIBMVolCons( self._params["softID"], self._params["kappaV"])) + def current_volume(self): + """ + Query the current volume of the soft object associated to this bond. + The volume is initialized once all :class:`IBM_Triel` bonds have + been added and the forces have been recalculated. + """ + return immersed_boundaries.get_current_volume(self._params['softID']) + class OifGlobalForces(BondedInteraction): diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index cd15ecc4409..a20c73ee4bb 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -142,6 +142,7 @@ python_test(FILE domain_decomposition.py MAX_NUM_PROC 4) python_test(FILE integrator_npt.py MAX_NUM_PROC 4) python_test(FILE integrator_npt_stats.py MAX_NUM_PROC 4 LABELS long) python_test(FILE integrator_steepest_descent.py MAX_NUM_PROC 4) +python_test(FILE ibm.py MAX_NUM_PROC 2) python_test(FILE dipolar_mdlc_p3m_scafacos_p2nfft.py MAX_NUM_PROC 1) python_test(FILE dipolar_direct_summation.py MAX_NUM_PROC 1 LABELS gpu) python_test(FILE dipolar_p3m.py MAX_NUM_PROC 2) diff --git a/testsuite/python/ibm.py b/testsuite/python/ibm.py new file mode 100644 index 00000000000..6eefe0ee303 --- /dev/null +++ b/testsuite/python/ibm.py @@ -0,0 +1,246 @@ +# +# Copyright (C) 2013-2019 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +import numpy as np +import itertools +import unittest as ut +import unittest_decorators as utx + +import espressomd +import espressomd.lb +import espressomd.interactions +try: + from espressomd.virtual_sites import VirtualSitesInertialessTracers +except ImportError: + pass + + +@utx.skipIfMissingFeatures(['VIRTUAL_SITES_INERTIALESS_TRACERS']) +class IBM(ut.TestCase): + '''Test IBM implementation with a Langevin thermostat.''' + system = espressomd.System(box_l=3 * [8.]) + system.time_step = 0.06 + system.cell_system.skin = 0.1 + + def tearDown(self): + self.system.part.clear() + self.system.actors.clear() + self.system.thermostat.turn_off() + + def compute_dihedral_angle(self, pos0, pos1, pos2, pos3): + # first normal vector + n1 = np.cross((pos1 - pos0), (pos2 - pos0)) + n2 = np.cross((pos2 - pos0), (pos3 - pos0)) + + norm1 = np.linalg.norm(n1) + norm2 = np.linalg.norm(n2) + n1 = n1 / norm1 + n2 = n2 / norm2 + + cos_alpha = min(1, np.dot(n1, n2)) + alpha = np.arccos(cos_alpha) + return alpha + + def test_tribend(self): + # two triangles with bending interaction + # move nodes, should relax back + + system = self.system + system.virtual_sites = VirtualSitesInertialessTracers() + system.thermostat.set_langevin(kT=0, gamma=10, seed=1) + + # Add four particles + p0 = system.part.add(pos=[4, 4, 4]) + p1 = system.part.add(pos=[4, 4, 5]) + p2 = system.part.add(pos=[4, 5, 5]) + p3 = system.part.add(pos=[4, 5, 4]) + + # Add first triel, weak modulus + tri1 = espressomd.interactions.IBM_Triel( + ind1=p0.id, ind2=p1.id, ind3=p2.id, elasticLaw="Skalak", k1=0.1, k2=0, maxDist=2.4) + system.bonded_inter.add(tri1) + p0.add_bond((tri1, p1, p2)) + + # Add second triel, strong modulus + tri2 = espressomd.interactions.IBM_Triel( + ind1=p0.id, ind2=p2.id, ind3=p3.id, elasticLaw="Skalak", k1=10, k2=0, maxDist=2.4) + system.bonded_inter.add(tri2) + p0.add_bond((tri2, p2, p3)) + + # Add bending + tribend = espressomd.interactions.IBM_Tribend( + ind1=p0.id, ind2=p1.id, ind3=p2.id, ind4=p3.id, kb=1, refShape="Initial") + system.bonded_inter.add(tribend) + p0.add_bond((tribend, p1, p2, p3)) + + # twist + system.part[:].pos = system.part[:].pos + np.random.random((4, 3)) + + # Perform integration + system.integrator.run(200) + angle = self.compute_dihedral_angle(p0.pos, p1.pos, p2.pos, p3.pos) + self.assertLess(angle, 2E-2) + + # IBM doesn't implement energy and pressure kernels. + energy = self.system.analysis.energy() + pressure = self.system.analysis.pressure() + self.assertAlmostEqual(energy['bonded'], 0., delta=1e-10) + self.assertAlmostEqual(pressure['bonded'], 0., delta=1e-10) + + def test_triel(self): + system = self.system + system.virtual_sites = VirtualSitesInertialessTracers() + system.thermostat.set_langevin(kT=0, gamma=1, seed=1) + + # Add particles: 0-2 are not bonded, 3-5 are bonded + non_bound = system.part.add(pos=[[4, 4, 4], [4, 4, 5], [4, 5, 5]]) + + p3 = system.part.add(pos=[1, 4, 4]) + p4 = system.part.add(pos=[1, 4, 5]) + p5 = system.part.add(pos=[1, 5, 5]) + + # Add triel for 3-5 + tri = espressomd.interactions.IBM_Triel( + ind1=p3.id, ind2=p4.id, ind3=p5.id, elasticLaw="Skalak", k1=15, + k2=0, maxDist=2.4) + system.bonded_inter.add(tri) + p3.add_bond((tri, p4, p5)) + + system.part[:].pos = system.part[:].pos + np.array(( + (0, 0, 0), (1, -.2, .3), (1, 1, 1), + (0, 0, 0), (1, -.2, .3), (1, 1, 1))) + + distorted_pos = np.copy(non_bound.pos) + + system.integrator.run(110) + dist1bound = system.distance(p3, p4) + dist2bound = system.distance(p3, p5) + + # check bonded particles. Distance should restore to initial config + self.assertAlmostEqual(dist1bound, 1, delta=0.05) + self.assertAlmostEqual(dist2bound, np.sqrt(2), delta=0.05) + + # check not bonded particles. Positions should still be distorted + np.testing.assert_allclose(np.copy(non_bound.pos), distorted_pos) + + def test_volcons(self): + '''Check volume conservation forces on a simple mesh (cube).''' + system = self.system + system.virtual_sites = VirtualSitesInertialessTracers() + system.thermostat.set_langevin(kT=0, gamma=1, seed=1) + + # Place particles on a cube. + positions = list(itertools.product((0, 1), repeat=3)) + positions = positions[:4] + positions[6:] + positions[4:6] + positions = np.array(positions) - 0.5 + mesh_center_ref = np.copy(system.box_l) / 2. + system.part.add(pos=positions + mesh_center_ref, virtual=8 * [True]) + + # Divide the cube. All triangle normals must point inside the mesh. + # Use the right hand rule to determine the order of the indices. + triangles = [(0, 1, 2), (1, 3, 2), + (2, 3, 4), (3, 5, 4), + (4, 5, 6), (5, 7, 6), + (6, 7, 0), (7, 1, 0), + (0, 2, 4), (0, 4, 6), + (1, 5, 3), (1, 7, 5)] + # Add triangle bonds that don't contribute to the force (infinite + # elasticity). These bonds are needed to calculate the mesh volume. + for id1, id2, id3 in triangles: + bond = espressomd.interactions.IBM_Triel( + ind1=id3, ind2=id2, ind3=id1, elasticLaw="Skalak", k1=0., k2=0., maxDist=3) + system.bonded_inter.add(bond) + system.part[id1].add_bond((bond, id2, id3)) + + # Add volume conservation force. + volCons = espressomd.interactions.IBM_VolCons(softID=15, kappaV=1.) + system.bonded_inter.add(volCons) + for p in system.part: + p.add_bond((volCons,)) + + # Run the integrator to initialize the mesh reference volume. + system.integrator.run(0, recalc_forces=True) + self.assertAlmostEqual(volCons.current_volume(), 1., delta=1e-10) + + # The restorative force is zero at the moment. + np.testing.assert_almost_equal(np.copy(self.system.part[:].f), 0.) + + # Double the cube dimensions. The volume increases by a factor of 8. + system.part[:].pos = 2. * positions + mesh_center_ref + system.integrator.run(0, recalc_forces=True) + self.assertAlmostEqual(volCons.current_volume(), 8., delta=1e-10) + + # Reference forces for that particular mesh geometry. + ref_forces = 1.75 * np.array( + [(1, 2, 2), (2, 1, -2), (2, -1, 1), (1, -2, -1), + (-1, -2, 2), (-2, -1, -2), (-2, 1, 1), (-1, 2, -1)]) + np.testing.assert_almost_equal( + np.copy(self.system.part[:].f), ref_forces) + + # IBM doesn't implement energy and pressure kernels. + energy = self.system.analysis.energy() + pressure = self.system.analysis.pressure() + self.assertAlmostEqual(energy['bonded'], 0., delta=1e-10) + self.assertAlmostEqual(pressure['bonded'], 0., delta=1e-10) + + # Add unthermalized LB. + system.thermostat.turn_off() + lbf = espressomd.lb.LBFluid( + kT=0.0, agrid=2, dens=1, visc=1.8, tau=system.time_step) + system.actors.add(lbf) + system.thermostat.set_lb(LB_fluid=lbf, act_on_virtual=False, gamma=1) + + # Check the cube is shrinking. The geometrical center shouldn't move. + volume_diff_ref = 5.805e-4 + # warmup + system.integrator.run(80) + # sampling + previous_volume = volCons.current_volume() + for _ in range(10): + system.integrator.run(20) + current_volume = volCons.current_volume() + volume_diff = previous_volume - current_volume + self.assertLess(current_volume, previous_volume) + self.assertAlmostEqual(volume_diff, volume_diff_ref, places=6) + previous_volume = current_volume + mesh_center = np.mean(self.system.part[:].pos, axis=0) + np.testing.assert_allclose(mesh_center, mesh_center_ref, rtol=1e-3) + # Halve the cube dimensions. The volume decreases by a factor of 8. + system.part[:].pos = 0.5 * positions + mesh_center_ref + system.integrator.run(0, recalc_forces=True) + self.assertAlmostEqual(volCons.current_volume(), 1. / 8., delta=1e-10) + + # Check the cube is expanding. The geometrical center shouldn't move. + volume_diff_ref = 6.6066e-7 + # warmup + system.integrator.run(80) + # sampling + previous_volume = volCons.current_volume() + for _ in range(10): + system.integrator.run(20) + current_volume = volCons.current_volume() + volume_diff = current_volume - previous_volume + self.assertGreater(current_volume, previous_volume) + self.assertAlmostEqual(volume_diff, volume_diff_ref, places=7) + previous_volume = current_volume + mesh_center = np.mean(self.system.part[:].pos, axis=0) + np.testing.assert_allclose(mesh_center, mesh_center_ref, rtol=1e-3) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/virtual_sites_tracers_common.py b/testsuite/python/virtual_sites_tracers_common.py index fc0ed9f0817..9717a5ac114 100644 --- a/testsuite/python/virtual_sites_tracers_common.py +++ b/testsuite/python/virtual_sites_tracers_common.py @@ -16,11 +16,8 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . # -import numpy as np - import espressomd from espressomd import shapes, lbboundaries -import espressomd.interactions try: from espressomd.virtual_sites import VirtualSitesInertialessTracers, VirtualSitesOff except ImportError: @@ -92,96 +89,6 @@ def test_advection(self): dist = self.lbf.get_interpolated_velocity(p.pos)[0] * system.time self.assertAlmostEqual(p.pos[0] / dist, 1, delta=0.005) - def compute_dihedral_angle(self, pos0, pos1, pos2, pos3): - # first normal vector - n1 = np.cross((pos1 - pos0), (pos2 - pos0)) - n2 = np.cross((pos2 - pos0), (pos3 - pos0)) - - norm1 = np.linalg.norm(n1) - norm2 = np.linalg.norm(n2) - n1 = n1 / norm1 - n2 = n2 / norm2 - - cos_alpha = min(1, np.dot(n1, n2)) - alpha = np.arccos(cos_alpha) - return alpha - - def test_tribend(self): - # two triangles with bending interaction - # move nodes, should relax back - - system = self.system - system.virtual_sites = VirtualSitesInertialessTracers() - system.thermostat.set_langevin(kT=0, gamma=10, seed=1) - - # Add four particles - p0 = system.part.add(pos=[5, 5, 5]) - p1 = system.part.add(pos=[5, 5, 6]) - p2 = system.part.add(pos=[5, 6, 6]) - p3 = system.part.add(pos=[5, 6, 5]) - - # Add first triel, weak modulus - tri1 = espressomd.interactions.IBM_Triel( - ind1=p0.id, ind2=p1.id, ind3=p2.id, elasticLaw="Skalak", k1=0.1, k2=0, maxDist=2.4) - system.bonded_inter.add(tri1) - p0.add_bond((tri1, p1, p2)) - - # Add second triel, strong modulus - tri2 = espressomd.interactions.IBM_Triel( - ind1=p0.id, ind2=p2.id, ind3=p3.id, elasticLaw="Skalak", k1=10, k2=0, maxDist=2.4) - system.bonded_inter.add(tri2) - p0.add_bond((tri2, p2, p3)) - - # Add bending - tribend = espressomd.interactions.IBM_Tribend( - ind1=p0.id, ind2=p1.id, ind3=p2.id, ind4=p3.id, kb=1, refShape="Initial") - system.bonded_inter.add(tribend) - p0.add_bond((tribend, p1, p2, p3)) - - # twist - system.part[:].pos = system.part[:].pos + np.random.random((4, 3)) - - # Perform integration - system.integrator.run(200) - angle = self.compute_dihedral_angle(p0.pos, p1.pos, p2.pos, p3.pos) - self.assertLess(angle, 2E-2) - - def test_triel(self): - system = self.system - system.virtual_sites = VirtualSitesInertialessTracers() - system.thermostat.set_langevin(kT=0, gamma=1, seed=1) - - # Add particles: 0-2 are not bonded, 3-5 are bonded - non_bound = system.part.add(pos=[[5, 5, 5], [5, 5, 6], [5, 6, 6]]) - - p3 = system.part.add(pos=[2, 5, 5]) - p4 = system.part.add(pos=[2, 5, 6]) - p5 = system.part.add(pos=[2, 6, 6]) - - # Add triel for 3-5 - tri = espressomd.interactions.IBM_Triel( - ind1=p3.id, ind2=p4.id, ind3=p5.id, elasticLaw="Skalak", k1=15, - k2=0, maxDist=2.4) - system.bonded_inter.add(tri) - p3.add_bond((tri, p4, p5)) - - system.part[:].pos = system.part[:].pos + np.array(( - (0, 0, 0), (1, -.2, .3), (1, 1, 1), - (0, 0, 0), (1, -.2, .3), (1, 1, 1))) - - distorted_pos = np.copy(non_bound.pos) - - system.integrator.run(110) - dist1bound = system.distance(p3, p4) - dist2bound = system.distance(p3, p5) - - # check bonded particles. Distance should restore to initial config - self.assertAlmostEqual(dist1bound, 1, delta=0.05) - self.assertAlmostEqual(dist2bound, np.sqrt(2), delta=0.05) - - # check not bonded particles. Positions should still be distorted - np.testing.assert_allclose(np.copy(non_bound.pos), distorted_pos) - def test_zz_without_lb(self): """Check behaviour without lb. Ignore non-virtual particles, complain on virtual ones.