Skip to content

Commit

Permalink
Rename Stress Tensor to Pressure Tensor
Browse files Browse the repository at this point in the history
  • Loading branch information
jngrad committed May 21, 2020
1 parent 969d976 commit c5a19f4
Show file tree
Hide file tree
Showing 41 changed files with 373 additions and 355 deletions.
34 changes: 17 additions & 17 deletions doc/sphinx/analysis.rst
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ The direct analysis commands can be classified into two types:
- :ref:`Center of mass`
- :ref:`Moment of inertia matrix`
- :ref:`Gyration tensor`
- :ref:`Stress Tensor`
- :ref:`Pressure Tensor`

.. _Energies:

Expand All @@ -41,7 +41,7 @@ Energies
:meth:`espressomd.analyze.Analysis.energy`

Returns the energies of the system.
The different energetic contributions to the total energy can also be obtained (kinetic, bonded,non-bonded, Coulomb).
The different energetic contributions to the total energy can also be obtained (kinetic, bonded, non-bonded, Coulomb).

For example, ::

Expand Down Expand Up @@ -203,26 +203,26 @@ constraints of any kind are not currently accounted for in the pressure
calculations. The pressure is no longer correct, e.g., when particles
are confined to a plane.

Note: The different contributions which are returned are the summands that arise from force splitting :math:`\vec{F}_{i,j}={\vec{F}_{i,j}}_\text{bonded}+{\vec{F}_{i,j}}_\text{nonbonded}+...` in the virial pressure formula. Later when the user calculates the ensemble average via e.g. :math:`\langle p \rangle \approx 1/N \sum_{i=1}^N p_i` however the ensemble average with all interactions present is performed. That means the contributions are not easy to interpret! Those are the contributions to the stress/pressure in a system where all interactions are present and therefore in a coupled system.
Note: The different contributions which are returned are the summands that arise from force splitting :math:`\vec{F}_{i,j}={\vec{F}_{i,j}}_\text{bonded}+{\vec{F}_{i,j}}_\text{nonbonded}+...` in the virial pressure formula. Later when the user calculates the ensemble average via e.g. :math:`\langle p \rangle \approx 1/N \sum_{i=1}^N p_i` however the ensemble average with all interactions present is performed. That means the contributions are not easy to interpret! Those are the contributions to the pressure in a system where all interactions are present and therefore in a coupled system.

.. _Stress Tensor:
.. _Pressure Tensor:

Stress Tensor
~~~~~~~~~~~~~
:meth:`espressomd.analyze.Analysis.stress_tensor`
Pressure Tensor
~~~~~~~~~~~~~~~
:meth:`espressomd.analyze.Analysis.pressure_tensor`

Computes the volume averaged instantaneous stress tensor of the system with options which are
described by in :meth:`espressomd.analyze.Analysis.stress_tensor`. It is called a stress tensor but the sign convention follows that of a pressure tensor.
In general do only use it for (on average) homogeneous systems. For inhomogeneous systems you need to use the local stress tensor.
Computes the volume averaged instantaneous pressure tensor of the system with options which are
described by in :meth:`espressomd.analyze.Analysis.pressure_tensor`.
In general do only use it for (on average) homogeneous systems. For inhomogeneous systems you need to use the local pressure tensor.

The instantaneous virial stress tensor is calculated by
The instantaneous virial pressure tensor is calculated by

.. math:: p_{(k,l)} = \frac{\sum_{i} {m_{i}v_{i}^{(k)}v_{i}^{(l)}}}{V} + \frac{\sum_{j>i}{F_{ij}^{(k)}r_{ij}^{(l)}}}{V}

where the notation is the same as for the pressure. The superscripts :math:`k`
and :math:`l` correspond to the components in the tensors and vectors.

If electrostatic interactions are present then also the coulombic parts of the stress tensor need to be calculated. If P3M is present, then the instantaneous stress tensor is added to the above equation in accordance with :cite:`essmann95a` :
If electrostatic interactions are present then also the coulombic parts of the pressure tensor need to be calculated. If P3M is present, then the instantaneous pressure tensor is added to the above equation in accordance with :cite:`essmann95a` :

.. math :: p^\text{Coulomb, P3M}_{(k,l)} =p^\text{Coulomb, P3M, dir}_{(k,l)} + p^\text{Coulomb, P3M, rec}_{(k,l)},
Expand All @@ -239,12 +239,12 @@ The long ranged (k-space) part is given by:
where :math:`S(\vec{k})` is the Fourier transformed charge density. Compared to Essmann we do not have the contribution :math:`p^\text{corr}_{k,l}` since we want to calculate the pressure that arises from all particles in the system.

Note: The different contributions which are returned are the summands that arise from force splitting :math:`\vec{F}_{i,j}={\vec{F}_{i,j}}_\text{bonded}+{\vec{F}_{i,j}}_\text{nonbonded}+...` in the virial stress tensor formula.
Later when the user calculates the stress tensor via :math:`\langle p_{(k,l)}\rangle \approx 1/N \sum_{i=1}^N p_{k,l}` however the ensemble average with all interactions present is performed.
That means the contributions are not easy to interpret! Those are the contributions to the stress/pressure in a system where all interactions are present and therefore in a coupled system.
Note: The different contributions which are returned are the summands that arise from force splitting :math:`\vec{F}_{i,j}={\vec{F}_{i,j}}_\text{bonded}+{\vec{F}_{i,j}}_\text{nonbonded}+...` in the virial pressure tensor formula.
Later when the user calculates the pressure tensor via :math:`\langle p_{(k,l)}\rangle \approx 1/N \sum_{i=1}^N p_{k,l}` however the ensemble average with all interactions present is performed.
That means the contributions are not easy to interpret! Those are the contributions to the pressure in a system where all interactions are present and therefore in a coupled system.

Note that the angular velocities of the particles are not included in
the calculation of the stress tensor.
the calculation of the pressure tensor.

.. _Chains:

Expand Down Expand Up @@ -504,7 +504,7 @@ documentation for all available observables in :mod:`espressomd.observables`.

- :class:`~espressomd.observables.Pressure`: Total scalar pressure (see :ref:`Pressure`)

- :class:`~espressomd.observables.PressureTensor`: Total pressure tensor (see :ref:`Stress Tensor`)
- :class:`~espressomd.observables.PressureTensor`: Total pressure tensor (see :ref:`Pressure Tensor`)

- :class:`~espressomd.observables.DPDStress`

Expand Down
12 changes: 6 additions & 6 deletions doc/sphinx/lb.rst
Original file line number Diff line number Diff line change
Expand Up @@ -181,12 +181,12 @@ Reading and setting properties of single lattice nodes

Appending three indices to the ``lb`` object returns an object that represents the selected LB grid node and allows one to access all of its properties::

lb[x, y, z].density # fluid density (one scalar for LB and CUDA)
lb[x, y, z].velocity # fluid velocity (a numpy array of three floats)
lb[x, y, z].stress # fluid pressure tensor (a symmetric 3x3 numpy array of floats)
lb[x, y, z].stress_neq # nonequilbrium part of the pressure tensor (as above)
lb[x, y, z].boundary # flag indicating whether the node is fluid or boundary (fluid: boundary=0, boundary: boundary != 0)
lb[x, y, z].population # 19 LB populations (a numpy array of 19 floats, check order from the source code)
lb[x, y, z].density # fluid density (one scalar for LB and CUDA)
lb[x, y, z].velocity # fluid velocity (a numpy array of three floats)
lb[x, y, z].pressure_tensor # fluid pressure tensor (a symmetric 3x3 numpy array of floats)
lb[x, y, z].pressure_tensor_neq # nonequilibrium part of the pressure tensor (as above)
lb[x, y, z].boundary # flag indicating whether the node is fluid or boundary (fluid: boundary=0, boundary: boundary != 0)
lb[x, y, z].population # 19 LB populations (a numpy array of 19 floats, check order from the source code)

All of these properties can be read and used in further calculations. Only the property ``population`` can be modified. The indices ``x,y,z`` are integers and enumerate the LB nodes in the three directions, starts with 0. To modify ``boundary``, refer to :ref:`Setting up boundary conditions`.

Expand Down
2 changes: 1 addition & 1 deletion doc/sphinx/particles.rst
Original file line number Diff line number Diff line change
Expand Up @@ -398,7 +398,7 @@ Please note:
needs to be adapted.

- The presence of rigid bodies constructed by means of virtual sites
adds a contribution to the pressure and stress tensor.
adds a contribution to the scalar pressure and pressure tensor.

.. _Inertialess lattice-Boltzmann tracers:

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -300,8 +300,8 @@
"The possible properties are:\n",
"\n",
"+ <tt>velocity</tt>: the fluid velocity (list of three floats)\n",
"+ <tt>stress</tt>: the pressure tensor (list of six floats: $\\Pi_{xx}$, $\\Pi_{xy}$, $\\Pi_{yy}$, $\\Pi_{xz}$, $\\Pi_{yz}$, and $\\Pi_{zz}$)\n",
"+ <tt>stress_neq</tt>: the nonequilibrium part of the pressure tensor, components as above.\n",
"+ <tt>pressure_tensor</tt>: the pressure tensor (3x3 matrix)\n",
"+ <tt>pressure_tensor_neq</tt>: the nonequilibrium part of the pressure tensor (3x3 matrix).\n",
"+ <tt>population</tt>: the 19 populations of the D3Q19 lattice.\n",
"+ <tt>boundary</tt>: the boundary flag.\n",
"+ <tt>density</tt>: the local density."
Expand Down
4 changes: 2 additions & 2 deletions src/core/Observable_stat.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -149,8 +149,8 @@ class Observable_stat : public Observable_stat_base {
}
};

/** Structure used only in the pressure and stress tensor calculation to
* distinguish non-bonded intra- and inter- molecular contributions.
/** Structure used only in the pressure calculation to distinguish
* non-bonded intra- and inter- molecular contributions.
*/
class Observable_stat_non_bonded : public Observable_stat_base {
public:
Expand Down
4 changes: 2 additions & 2 deletions src/core/bonded_interactions/bonded_interactions.dox
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
* - write functions for calculating force and energy
* - add calls to the force and energy calculation functions to the force
* calculation in the integration loop as well as to energy and
* pressure/stress tensor analysis
* pressure tensor analysis
* * Python interface:
* - import the definition of the interaction struct from the core
* - implement a class for the bonded interaction derived from the Python
Expand Down Expand Up @@ -146,7 +146,7 @@
* retval = fene_pair_energy(iaparams, dx);
* break;
* @endcode
* * Pressure, stress tensor and virial calculation: if your bonded
* * Pressure tensor and virial calculation: if your bonded
* interaction is not a pair bond or modifies the particles involved,
* you have to implement a custom solution for virial calculation.
* The pressure occurs twice, once for the parallelized isotropic pressure
Expand Down
2 changes: 1 addition & 1 deletion src/core/communication.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -220,7 +220,7 @@ void mpi_bcast_max_seen_particle_type(int s);
* \arg for \ref GatherStats::energy, calculate and reduce (sum up)
* energies, using \ref energy_calc.
* \arg for \ref GatherStats::pressure, calculate and reduce (sum up)
* pressure, stress tensor, using \ref pressure_calc.
* pressure, using \ref pressure_calc.
* \arg for \ref GatherStats::pressure_v_comp, calculate and reduce
* (sum up) instantaneous pressure, using \ref pressure_calc.
* \arg for \ref GatherStats::lb_fluid_momentum, use
Expand Down
6 changes: 3 additions & 3 deletions src/core/electrostatics_magnetostatics/coulomb.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,9 @@ void calc_pressure_long_range(Observable_stat &virials,
break;
case COULOMB_P3M: {
p3m_charge_assign(particles);
auto const p3m_stress = p3m_calc_kspace_stress();
std::copy_n(p3m_stress.data(), 9, p_tensor.coulomb.begin() + 9);
virials.coulomb[1] = p3m_stress[0] + p3m_stress[4] + p3m_stress[8];
auto const p3m_p_tensor = p3m_calc_kspace_pressure_tensor();
std::copy_n(p3m_p_tensor.data(), 9, p_tensor.coulomb.begin() + 9);
virials.coulomb[1] = p3m_p_tensor[0] + p3m_p_tensor[4] + p3m_p_tensor[8];

break;
}
Expand Down
32 changes: 16 additions & 16 deletions src/core/electrostatics_magnetostatics/p3m.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -419,13 +419,13 @@ double dipole_correction_energy(Utils::Vector3d const &box_dipole) {
}
} // namespace

/** @details Calculate the long range electrostatics part of the stress
/** @details Calculate the long range electrostatics part of the pressure
* tensor. This is part \f$\Pi_{\textrm{dir}, \alpha, \beta}\f$ eq. (2.6)
* in @cite essmann95a. The part \f$\Pi_{\textrm{corr}, \alpha, \beta}\f$
* eq. (2.8) is not present here since M is the empty set in our simulations.
*/
Utils::Vector9d p3m_calc_kspace_stress() {
Utils::Vector9d node_k_space_stress{};
Utils::Vector9d p3m_calc_kspace_pressure_tensor() {
Utils::Vector9d node_k_space_pressure_tensor{};

if (p3m.sum_q2 > 0) {
p3m.sm.gather_grid(p3m.rs_mesh.data(), comm_cart, p3m.local_mesh.dim);
Expand Down Expand Up @@ -461,25 +461,25 @@ Utils::Vector9d p3m_calc_kspace_stress() {

diagonal += node_k_space_energy;
auto const prefactor = node_k_space_energy * vterm;
node_k_space_stress[0] += prefactor * kx * kx; /* sigma_xx */
node_k_space_stress[1] += prefactor * kx * ky; /* sigma_xy */
node_k_space_stress[2] += prefactor * kx * kz; /* sigma_xz */
node_k_space_stress[3] += prefactor * ky * kx; /* sigma_yx */
node_k_space_stress[4] += prefactor * ky * ky; /* sigma_yy */
node_k_space_stress[5] += prefactor * ky * kz; /* sigma_yz */
node_k_space_stress[6] += prefactor * kz * kx; /* sigma_zx */
node_k_space_stress[7] += prefactor * kz * ky; /* sigma_zy */
node_k_space_stress[8] += prefactor * kz * kz; /* sigma_zz */
node_k_space_pressure_tensor[0] += prefactor * kx * kx; /* sigma_xx */
node_k_space_pressure_tensor[1] += prefactor * kx * ky; /* sigma_xy */
node_k_space_pressure_tensor[2] += prefactor * kx * kz; /* sigma_xz */
node_k_space_pressure_tensor[3] += prefactor * ky * kx; /* sigma_yx */
node_k_space_pressure_tensor[4] += prefactor * ky * ky; /* sigma_yy */
node_k_space_pressure_tensor[5] += prefactor * ky * kz; /* sigma_yz */
node_k_space_pressure_tensor[6] += prefactor * kz * kx; /* sigma_zx */
node_k_space_pressure_tensor[7] += prefactor * kz * ky; /* sigma_zy */
node_k_space_pressure_tensor[8] += prefactor * kz * kz; /* sigma_zz */
}
}
}
node_k_space_stress[0] += diagonal;
node_k_space_stress[4] += diagonal;
node_k_space_stress[8] += diagonal;
node_k_space_pressure_tensor[0] += diagonal;
node_k_space_pressure_tensor[4] += diagonal;
node_k_space_pressure_tensor[8] += diagonal;
}

auto const force_prefac = coulomb.prefactor / (2.0 * box_geo.volume());
return force_prefac * node_k_space_stress;
return force_prefac * node_k_space_pressure_tensor;
}

double p3m_calc_kspace_forces(bool force_flag, bool energy_flag,
Expand Down
4 changes: 2 additions & 2 deletions src/core/electrostatics_magnetostatics/p3m.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -152,8 +152,8 @@ void p3m_scaleby_box_l();
double p3m_calc_kspace_forces(bool force_flag, bool energy_flag,
const ParticleRange &particles);

/** Compute the k-space part of the stress tensor **/
Utils::Vector9d p3m_calc_kspace_stress();
/** Compute the k-space part of the pressure tensor **/
Utils::Vector9d p3m_calc_kspace_pressure_tensor();

/** Sanity checks */
bool p3m_sanity_checks();
Expand Down
6 changes: 3 additions & 3 deletions src/core/grid_based_algorithms/lb.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1279,9 +1279,9 @@ Utils::Vector3d lb_calc_momentum_density(std::array<double, 19> const &modes,
modes[3] + 0.5 * force_density[2]}};
}

Utils::Vector6d lb_calc_stress(std::array<double, 19> const &modes,
Utils::Vector3d const &force_density,
const LB_Parameters &lb_parameters) {
Utils::Vector6d lb_calc_pressure_tensor(std::array<double, 19> const &modes,
Utils::Vector3d const &force_density,
const LB_Parameters &lb_parameters) {
auto const momentum_density = lb_calc_momentum_density(modes, force_density);
auto const density = lb_calc_density(modes, lb_parameters);
using Utils::sqr;
Expand Down
10 changes: 5 additions & 5 deletions src/core/grid_based_algorithms/lb.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,8 @@
* This increases mainly cache efficiency. This is achieved by
* @ref lb_collide_stream.
*
* The hydrodynamic fields, corresponding to density, velocity and stress, are
* stored in @ref LB_FluidNode in the array @ref lbfields, the populations
* The hydrodynamic fields, corresponding to density, velocity and pressure,
* are stored in @ref LB_FluidNode in the array @ref lbfields, the populations
* in @ref LB_Fluid in the array @ref lbfluid which is constructed as
* 2 x (Nx x Ny x Nz) x 19 array.
*
Expand Down Expand Up @@ -206,9 +206,9 @@ double lb_calc_density(std::array<double, 19> const &modes,
const LB_Parameters &lb_parameters);
Utils::Vector3d lb_calc_momentum_density(std::array<double, 19> const &modes,
Utils::Vector3d const &force_density);
Utils::Vector6d lb_calc_stress(std::array<double, 19> const &modes,
Utils::Vector3d const &force_density,
const LB_Parameters &lb_parameters);
Utils::Vector6d lb_calc_pressure_tensor(std::array<double, 19> const &modes,
Utils::Vector3d const &force_density,
const LB_Parameters &lb_parameters);

/** Calculation of hydrodynamic modes.
*
Expand Down
6 changes: 3 additions & 3 deletions src/core/grid_based_algorithms/lb_collective_interface.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -139,14 +139,14 @@ auto mpi_lb_get_momentum_density(Utils::Vector3i const &index) {

REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_momentum_density)

auto mpi_lb_get_stress(Utils::Vector3i const &index) {
auto mpi_lb_get_pressure_tensor(Utils::Vector3i const &index) {
return detail::lb_calc_fluid_kernel(
index, [&](auto modes, auto force_density) {
return lb_calc_stress(modes, force_density, lbpar);
return lb_calc_pressure_tensor(modes, force_density, lbpar);
});
}

REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_stress)
REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_pressure_tensor)

void mpi_bcast_lb_params_slave(LBParam field, LB_Parameters const &params) {
lbpar = params;
Expand Down
2 changes: 1 addition & 1 deletion src/core/grid_based_algorithms/lb_collective_interface.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ boost::optional<int> mpi_lb_get_boundary_flag(Utils::Vector3i const &index);
boost::optional<Utils::Vector3d>
mpi_lb_get_momentum_density(Utils::Vector3i const &index);
boost::optional<Utils::Vector6d>
mpi_lb_get_stress(Utils::Vector3i const &index);
mpi_lb_get_pressure_tensor(Utils::Vector3i const &index);

/* collective setter functions */
void mpi_lb_set_population(Utils::Vector3i const &index,
Expand Down
Loading

0 comments on commit c5a19f4

Please sign in to comment.