Skip to content

Commit

Permalink
LB: per particle gamma
Browse files Browse the repository at this point in the history
  • Loading branch information
@christophlohrmann @jngrad
christophlohrmann authored and jngrad committed Jun 13, 2023
1 parent 452a6e6 commit 3c2cb0b
Show file tree
Hide file tree
Showing 7 changed files with 286 additions and 132 deletions.
50 changes: 46 additions & 4 deletions src/core/grid_based_algorithms/lb_particle_coupling.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include "errorhandling.hpp"
#include "grid.hpp"
#include "random.hpp"
#include "thermostat.hpp"

#include "grid_based_algorithms/lb_interface.hpp"
#include "grid_based_algorithms/lb_interpolation.hpp"
Expand All @@ -35,6 +36,7 @@

#include <boost/mpi.hpp>

#include <cmath>
#include <cstdint>
#include <initializer_list>
#include <limits>
Expand Down Expand Up @@ -121,6 +123,21 @@ Utils::Vector3d lb_particle_coupling_drift_vel_offset(const Particle &p) {
return vel_offset;
}

static Thermostat::GammaType lb_handle_particle_anisotropy(Particle const &p) {
auto const lb_gamma = lb_lbcoupling_get_gamma();
#ifdef THERMOSTAT_PER_PARTICLE
auto const &partcl_gamma = p.gamma();
#ifdef PARTICLE_ANISOTROPY
auto const default_gamma = Thermostat::GammaType::broadcast(lb_gamma);
#else
auto const default_gamma = lb_gamma;
#endif // PARTICLE_ANISOTROPY
return handle_particle_gamma(partcl_gamma, default_gamma);
#else
return lb_gamma;
#endif // THERMOSTAT_PER_PARTICLE
}

Utils::Vector3d lb_drag_force(Particle const &p,
Utils::Vector3d const &shifted_pos,
Utils::Vector3d const &vel_offset) {
Expand All @@ -131,8 +148,10 @@ Utils::Vector3d lb_drag_force(Particle const &p,
LB::get_lattice_speed();

Utils::Vector3d v_drift = interpolated_u + vel_offset;
auto const gamma = lb_handle_particle_anisotropy(p);

/* calculate viscous force (eq. (9) @cite ahlrichs99a) */
return -lb_lbcoupling_get_gamma() * (p.v() - v_drift);
return Utils::hadamard_product(gamma, v_drift - p.v());
}

/**
Expand Down Expand Up @@ -210,7 +229,7 @@ void add_swimmer_force(Particle const &p, double time_step) {

namespace LB {

Utils::Vector3d ParticleCoupling::get_noise_term(int pid) const {
Utils::Vector3d ParticleCoupling::get_noise_term(Particle const &p) const {
if (not m_thermalized) {
return Utils::Vector3d{};
}
Expand All @@ -219,8 +238,14 @@ Utils::Vector3d ParticleCoupling::get_noise_term(int pid) const {
throw std::runtime_error(
"Access to uninitialized LB particle coupling RNG counter");
}
using std::sqrt;
using Utils::sqrt;
auto const counter = rng_counter->value();
return m_noise * Random::noise_uniform<RNGSalt::PARTICLES>(counter, 0, pid);
auto const gamma = lb_handle_particle_anisotropy(p);
return m_noise_pref_wo_gamma *
Utils::hadamard_product(
sqrt(gamma),
Random::noise_uniform<RNGSalt::PARTICLES>(counter, 0, p.id()));
}

void ParticleCoupling::kernel(Particle &p) {
Expand All @@ -233,7 +258,7 @@ void ParticleCoupling::kernel(Particle &p) {
if (in_local_halo(pos)) {
auto const vel_offset = lb_particle_coupling_drift_vel_offset(p);
auto const drag_force = lb_drag_force(p, pos, vel_offset);
auto const random_force = get_noise_term(p.id());
auto const random_force = get_noise_term(p);
coupling_force = drag_force + random_force;
break;
}
Expand All @@ -259,6 +284,20 @@ bool CouplingBookkeeping::is_ghost_for_local_particle(Particle const &p) const {
return not ::cell_structure.get_local_particle(p.id())->is_ghost();
}

#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY)
static void lb_coupling_sanity_checks(Particle const &p) {
/*
lb does (at the moment) not support rotational particle coupling.
Consequently, anisotropic particles are also not supported.
*/
auto const &p_gamma = p.gamma();
if (p_gamma[0] != p_gamma[1] or p_gamma[1] != p_gamma[2]) {
runtimeErrorMsg() << "anisotropic particle (id " << p.id()
<< ") coupled to LB.";
}
}
#endif

void couple_particles(bool couple_virtual, ParticleRange const &real_particles,
ParticleRange const &ghost_particles, double time_step) {
ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
Expand All @@ -269,6 +308,9 @@ void couple_particles(bool couple_virtual, ParticleRange const &real_particles,
for (auto const &particle_range : {real_particles, ghost_particles}) {
for (auto &p : particle_range) {
if (bookkeeping.should_be_coupled(p)) {
#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY)
lb_coupling_sanity_checks(p);
#endif
coupling.kernel(p);
}
}
Expand Down
22 changes: 9 additions & 13 deletions src/core/grid_based_algorithms/lb_particle_coupling.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,31 +139,27 @@ class ParticleCoupling {
bool m_couple_virtual;
bool m_thermalized;
double m_time_step;
double m_noise;
double m_noise_pref_wo_gamma;

public:
auto get_noise(double kT, double gamma) const {
/* Eq. (16) @cite ahlrichs99a.
* The factor 12 comes from the fact that we use random numbers
* from -0.5 to 0.5 (equally distributed) which have variance 1/12.
* The time step comes from the discretization.
*/
auto constexpr variance_correction = 12.;
return std::sqrt(variance_correction * 2. * gamma * kT / m_time_step);
}

ParticleCoupling(bool couple_virtual, double time_step, double kT)
: m_couple_virtual{couple_virtual}, m_thermalized{kT != 0.},
m_time_step{time_step} {
assert(kT >= 0.);
m_noise = get_noise(kT, lb_particle_coupling.gamma);
/* Eq. (16) @cite ahlrichs99a, without the gamma term.
* The factor 12 comes from the fact that we use random numbers
* from -0.5 to 0.5 (equally distributed) which have variance 1/12.
* The time step comes from the discretization.
*/
auto constexpr variance_inv = 12.;
m_noise_pref_wo_gamma = std::sqrt(variance_inv * 2. * kT / time_step);
}

ParticleCoupling(bool couple_virtual, double time_step)
: ParticleCoupling(couple_virtual, time_step,
LB::get_kT() * Utils::sqr(LB::get_lattice_speed())) {}

Utils::Vector3d get_noise_term(int pid) const;
Utils::Vector3d get_noise_term(Particle const &p) const;
void kernel(Particle &p);
};

Expand Down
65 changes: 50 additions & 15 deletions src/core/thermostat.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,9 @@
* Implementation in \ref thermostat.cpp.
*/

#include "Particle.hpp"
#include "rotation.hpp"

#include "config/config.hpp"

#include <utils/Counter.hpp>
Expand Down Expand Up @@ -55,16 +58,24 @@ using GammaType = double;
} // namespace Thermostat

namespace {
/** @name Integrators parameters sentinels.
* These functions return the sentinel value for the Langevin/Brownian
* parameters, indicating that they have not been set yet.
/**
* @brief Value for unset friction coefficient.
* Sentinel value for the Langevin/Brownian parameters,
* indicating that they have not been set yet.
*/
/**@{*/
constexpr double sentinel(double) { return -1.0; }
constexpr Utils::Vector3d sentinel(Utils::Vector3d) {
return {-1.0, -1.0, -1.0};
}
/**@}*/
#ifdef PARTICLE_ANISOTROPY
constexpr Thermostat::GammaType gamma_sentinel{{-1.0, -1.0, -1.0}};
#else
constexpr Thermostat::GammaType gamma_sentinel{-1.0};
#endif
/**
* @brief Value for a null friction coefficient.
*/
#ifdef PARTICLE_ANISOTROPY
constexpr Thermostat::GammaType gamma_null{{0.0, 0.0, 0.0}};
#else
constexpr Thermostat::GammaType gamma_null{0.0};
#endif
} // namespace

/************************************************
Expand All @@ -85,6 +96,30 @@ extern double temperature;
/** True if the thermostat should act on virtual particles. */
extern bool thermo_virtual;

#ifdef THERMOSTAT_PER_PARTICLE
inline auto const &
handle_particle_gamma(Thermostat::GammaType const &particle_gamma,
Thermostat::GammaType const &default_gamma) {
return particle_gamma >= gamma_null ? particle_gamma : default_gamma;
}
#endif

inline auto
handle_particle_anisotropy(Particle const &p,
Thermostat::GammaType const &gamma_body) {
#ifdef PARTICLE_ANISOTROPY
auto const aniso_flag =
(gamma_body[0] != gamma_body[1]) || (gamma_body[1] != gamma_body[2]);
const Utils::Matrix<double, 3, 3> gamma_matrix =
boost::qvm::diag_mat(gamma_body);
auto const gamma_space =
aniso_flag ? convert_body_to_space(p, gamma_matrix) : gamma_matrix;
return gamma_space;
#else
return gamma_body;
#endif
}

/************************************************
* parameter structs
************************************************/
Expand Down Expand Up @@ -142,9 +177,9 @@ struct LangevinThermostat : public BaseThermostat {
/** @name Parameters */
/**@{*/
/** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */
GammaType gamma = sentinel(GammaType{});
GammaType gamma = gamma_sentinel;
/** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */
GammaType gamma_rotation = sentinel(GammaType{});
GammaType gamma_rotation = gamma_sentinel;
/**@}*/
/** @name Prefactors */
/**@{*/
Expand Down Expand Up @@ -219,9 +254,9 @@ struct BrownianThermostat : public BaseThermostat {
/** @name Parameters */
/**@{*/
/** Translational friction coefficient @f$ \gamma_{\text{trans}} @f$. */
GammaType gamma = sentinel(GammaType{});
GammaType gamma = gamma_sentinel;
/** Rotational friction coefficient @f$ \gamma_{\text{rot}} @f$. */
GammaType gamma_rotation = sentinel(GammaType{});
GammaType gamma_rotation = gamma_sentinel;
/**@}*/
/** @name Prefactors */
/**@{*/
Expand All @@ -230,13 +265,13 @@ struct BrownianThermostat : public BaseThermostat {
* @f$ D_{\text{trans}} = k_B T/\gamma_{\text{trans}} @f$
* the translational diffusion coefficient.
*/
GammaType sigma_pos = sentinel(GammaType{});
GammaType sigma_pos = gamma_sentinel;
/** Rotational noise standard deviation.
* Stores @f$ \sqrt{2D_{\text{rot}}} @f$ with
* @f$ D_{\text{rot}} = k_B T/\gamma_{\text{rot}} @f$
* the rotational diffusion coefficient.
*/
GammaType sigma_pos_rotation = sentinel(GammaType{});
GammaType sigma_pos_rotation = gamma_sentinel;
/** Translational velocity noise standard deviation.
* Stores @f$ \sqrt{k_B T} @f$.
*/
Expand Down
56 changes: 16 additions & 40 deletions src/core/thermostats/langevin_inline.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,16 +50,13 @@ friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p,
}

// Determine prefactors for the friction and the noise term
Thermostat::GammaType pref_friction = langevin.pref_friction;
Thermostat::GammaType pref_noise = langevin.pref_noise;
#ifdef THERMOSTAT_PER_PARTICLE
// override default if particle-specific gamma
if (p.gamma() >= Thermostat::GammaType{}) {
auto const gamma =
p.gamma() >= Thermostat::GammaType{} ? p.gamma() : langevin.gamma;
pref_friction = -gamma;
pref_noise = LangevinThermostat::sigma(kT, time_step, gamma);
}
auto const gamma = handle_particle_gamma(p.gamma(), langevin.gamma);
auto const pref_friction = -gamma;
auto const pref_noise = LangevinThermostat::sigma(kT, time_step, gamma);
#else
auto const pref_friction = langevin.pref_friction;
auto const pref_noise = langevin.pref_noise;
#endif // THERMOSTAT_PER_PARTICLE

// Get effective velocity in the thermostatting
Expand All @@ -70,26 +67,9 @@ friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p,
#else
auto const &velocity = p.v();
#endif // ENGINE
#ifdef PARTICLE_ANISOTROPY
// Particle frictional isotropy check
auto const aniso_flag = (pref_friction[0] != pref_friction[1]) ||
(pref_friction[1] != pref_friction[2]);

// In case of anisotropic particle: body-fixed reference frame. Otherwise:
// lab-fixed reference frame.
const Utils::Matrix<double, 3, 3> fric_mat =
boost::qvm::diag_mat(pref_friction);
const Utils::Matrix<double, 3, 3> noise_mat =
boost::qvm::diag_mat(pref_noise);

auto const friction_op =
aniso_flag ? convert_body_to_space(p, fric_mat) : fric_mat;
auto const noise_op =
aniso_flag ? convert_body_to_space(p, noise_mat) : noise_mat;
#else
auto const &friction_op = pref_friction;
auto const &noise_op = pref_noise;
#endif // PARTICLE_ANISOTROPY
auto const friction_op = handle_particle_anisotropy(p, pref_friction);
auto const noise_op = handle_particle_anisotropy(p, pref_noise);

return friction_op * velocity +
noise_op * Random::noise_uniform<RNGSalt::LANGEVIN>(
Expand All @@ -111,23 +91,19 @@ friction_thermo_langevin_rotation(LangevinThermostat const &langevin,
Particle const &p, double time_step,
double kT) {

auto pref_friction = -langevin.gamma_rotation;
auto pref_noise = langevin.pref_noise_rotation;

#ifdef THERMOSTAT_PER_PARTICLE
// override default if particle-specific gamma
if (p.gamma_rot() >= Thermostat::GammaType{}) {
auto const gamma = p.gamma_rot() >= Thermostat::GammaType{}
? p.gamma_rot()
: langevin.gamma_rotation;
pref_friction = -gamma;
pref_noise = LangevinThermostat::sigma(kT, time_step, gamma);
}
auto const gamma =
handle_particle_gamma(p.gamma_rot(), langevin.gamma_rotation);
auto const pref_friction = gamma;
auto const pref_noise = LangevinThermostat::sigma(kT, time_step, gamma);
#else
auto const pref_friction = langevin.gamma_rotation;
auto const pref_noise = langevin.pref_noise_rotation;
#endif // THERMOSTAT_PER_PARTICLE

auto const noise = Random::noise_uniform<RNGSalt::LANGEVIN_ROT>(
langevin.rng_counter(), langevin.rng_seed(), p.id());
return hadamard_product(pref_friction, p.omega()) +
return -hadamard_product(pref_friction, p.omega()) +
hadamard_product(pref_noise, noise);
}

Expand Down
Loading

0 comments on commit 3c2cb0b

Please sign in to comment.