Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix Constraint checkpointing #3105

Merged
merged 9 commits into from
Aug 30, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 16 additions & 3 deletions src/core/field_coupling/fields/Interpolated.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -51,13 +51,20 @@ void deep_copy(boost::multi_array<T, 3> &dst,
} // namespace detail

/**
* @brief A vector field interpolated from a regular grid.
* @brief A vector or scalar field interpolated from a regular grid.
*
* This is an interpolation wrapper around a boost::multi_array,
* which can be evaluated on any point in space by spline interpolation.
* This is an interpolation wrapper around a boost::multi_array,
* which can be evaluated on any point in space by spline interpolation.
*
* @tparam T Underlying type of the field values, see @ref value_type
* @tparam codim Dimension of the field: 3 for a vector field,
* 1 for a scalar field.
*/
template <typename T, size_t codim> class Interpolated {
public:
/** Type of the values, usually @ref Utils::Vector<T, 3> for vector fields
* and @p T for scalar fields
*/
using value_type =
typename Utils::decay_to_scalar<Utils::Vector<T, codim>>::type;
using jacobian_type = detail::jacobian_type<T, codim>;
Expand Down Expand Up @@ -97,6 +104,12 @@ template <typename T, size_t codim> class Interpolated {
return {m_global_field.shape(), m_global_field.shape() + 3};
}

/** Serialize field */
std::vector<T> field_data_flat() const {
auto const *data = reinterpret_cast<T const *>(m_global_field.data());
return std::vector<T>(data, data + codim * m_global_field.num_elements());
}

/*
* @brief Evaluate f at pos with the field value as argument.
*/
Expand Down
144 changes: 95 additions & 49 deletions src/python/espressomd/constraints.py
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@ class ShapeBasedConstraint(Constraint):
only useful if penetrable is ``True``.
particle_type : int
Interaction type of the constraint.
particle_velocity : array of :obj:`float`
particle_velocity : array_like of :obj:`float`
Interaction velocity of the boundary
penetrable : :obj:`bool`
Whether particles are allowed to penetrate the constraint.
Expand Down Expand Up @@ -190,7 +190,7 @@ class HomogeneousMagneticField(Constraint):
"""
Attributes
----------
H : array of :obj:`float`
H : (3,) array_like of :obj:`float`
Magnetic field vector. Describes both field direction and
strength of the magnetic field (via length of the vector).

Expand All @@ -214,25 +214,38 @@ class _Interpolated(Constraint):
box, e.g. the most up right back point has to be at least at
box + 0.5 * grid_spacing. There are convenience functions on this
class that can calculate the required grid dimensions and the coordinates.
See also the examples on ForceField.

Arguments
----------
field : (M, N, O, P) array_like of :obj:`float`
The actual field on a grid of size (M, N, O) with dimension P.
grid_spacing : (3,) array_like of :obj:`float`
Spacing of the grid points.

Attributes
----------

field_data: array_like :obj:`float`
The actual field. Please be aware that depending on the interpolation
field : (M, N, O, P) array_like of :obj:`float`
The actual field on a grid of size (M, N, O) with dimension P.
Please be aware that depending on the interpolation
order additional points are used on the boundaries.

grid_spacing: array_like :obj:`float`
The spacing of the grid points.
grid_spacing : array_like of :obj:`float`
Spacing of the grid points.

"""
origin : (3,) array_like of :obj:`float`
Coordinates of the grid origin.

def __init__(self, field, **kwargs):
shape, codim = self._unpack_dims(field)
"""

super().__init__(_field_shape=shape, _field_codim=codim,
_field_data=field.flatten(), **kwargs)
def __init__(self, **kwargs):
if "oid" not in kwargs:
field = kwargs.pop("field")
shape, codim = self._unpack_dims(field)
super().__init__(_field_shape=shape, _field_codim=codim,
_field_data=field.flatten(), **kwargs)
else:
super().__init__(**kwargs)

@classmethod
def required_dims(cls, box_size, grid_spacing):
Expand Down Expand Up @@ -321,19 +334,23 @@ class ForceField(_Interpolated):
"""
A generic tabulated force field that applies a per-particle scaling factor.

Attributes
Arguments
----------
field : (M, N, O, 3) array_like of :obj:`float`
Forcefield amplitude on a grid of size (M, N, O).
grid_spacing : (3,) array_like of :obj:`float`
Spacing of the grid points.
default_scale : :obj:`float`
Scaling factor for particles that have no individual scaling factor.
particle_scales: array_like (:obj:`int`, :obj:`float`)
particle_scales : array_like of (:obj:`int`, :obj:`float`)
A list of tuples of ids and scaling factors. For
particles in the list the interaction is scaled with
their individual scaling factor before it is applied.

"""

def __init__(self, field, **kwargs):
super().__init__(field, **kwargs)
def __init__(self, **kwargs):
super().__init__(**kwargs)

_codim = 3
_so_name = "Constraints::ForceField"
Expand All @@ -343,25 +360,28 @@ def __init__(self, field, **kwargs):
class PotentialField(_Interpolated):

"""
A generic tabulated force field that applies a per particle
A generic tabulated force field that applies a per-particle
scaling factor. The forces are calculated numerically from
the data by finite differences. The potential is interpolated
from the provided data.

Attributes
Arguments
----------
field : (M, N, O, 1) array_like of :obj:`float`
Potential on a grid of size (M, N, O).
grid_spacing : (3,) array_like of :obj:`float`
Spacing of the grid points.
default_scale : :obj:`float`
Scaling factor for particles that have no
individual scaling factor.
particle_scales: array_like (:obj:`int`, :obj:`float`)
Scaling factor for particles that have no individual scaling factor.
particle_scales : array_like (:obj:`int`, :obj:`float`)
A list of tuples of ids and scaling factors. For
particles in the list the interaction is scaled with
their individual scaling factor before it is applied.

"""

def __init__(self, field, **kwargs):
super().__init__(field, **kwargs)
def __init__(self, **kwargs):
super().__init__(**kwargs)

_codim = 1
_so_name = "Constraints::PotentialField"
Expand All @@ -375,15 +395,17 @@ class Gravity(Constraint):

:math:`F = m \\cdot g`

Attributes
Arguments
----------
g : array of :obj:`float`
g : (3,) array_like of :obj:`float`
The gravitational acceleration.

"""

def __init__(self, g):
super().__init__(value=g)
def __init__(self, **kwargs):
if "oid" not in kwargs:
kwargs["value"] = kwargs.pop("g")
super().__init__(**kwargs)

@property
def g(self):
Expand All @@ -408,18 +430,21 @@ class LinearElectricPotential(Constraint):

where :math:`q` is the charge of the particle.

Attributes
Arguments
----------
E : array of :obj:`float`
E : array_like of :obj:`float`
The electric field.

phi0 : :obj:`float`
The potential at the origin

"""

def __init__(self, E, phi0=0):
super().__init__(A=-E, b=phi0)
def __init__(self, phi0=0, **kwargs):
if "oid" not in kwargs:
kwargs["A"] = -np.array(kwargs.pop("E"))
kwargs["b"] = phi0
super().__init__(**kwargs)

@property
def E(self):
Expand Down Expand Up @@ -448,24 +473,28 @@ class ElectricPlaneWave(Constraint):
This can be used to generate a homogeneous AC
field by setting k to zero.

Attributes
Arguments
----------
E0 : array of :obj:`float`
The amplitude of the electric field.
k : array of :obj:`float`
E0 : array_like of :obj:`float`
Amplitude of the electric field.
k : array_like of :obj:`float`
Wave vector of the wave
omega : :obj:`float`
Frequency of the wave
phi : :obj:`float`
Optional phase shift, defaults to 0.
phi : :obj:`float`, optional
Phase shift

"""

_so_name = "Constraints::ElectricPlaneWave"

def __init__(self, E0, k, omega, phi=0):
super().__init__(amplitude=E0, wave_vector=k, frequency=omega,
phase=phi)
def __init__(self, phi=0, **kwargs):
if "oid" not in kwargs:
kwargs["amplitude"] = kwargs.pop("E0")
kwargs["wave_vector"] = kwargs.pop("k")
kwargs["frequency"] = kwargs.pop("omega")
kwargs["phase"] = phi
super().__init__(**kwargs)

@property
def E0(self):
Expand Down Expand Up @@ -495,10 +524,19 @@ class FlowField(_Interpolated):

where :math:`v` is the velocity of the particle.

Arguments
----------
field : (M, N, O, 3) array_like of :obj:`float`
Field velocity on a grid of size (M, N, O)
grid_spacing : (3,) array_like of :obj:`float`
Spacing of the grid points.
gamma : :obj:`float`
Coupling constant

"""

def __init__(self, field, **kwargs):
super().__init__(field, **kwargs)
def __init__(self, **kwargs):
super().__init__(**kwargs)

_codim = 3
_so_name = "Constraints::FlowField"
Expand All @@ -517,15 +555,17 @@ class HomogeneousFlowField(Constraint):

Attributes
----------
gamma : :obj:`float`
The coupling constant
u : (3,) array_like of :obj:`float`
The velocity of the field.
Field velocity.
gamma : :obj:`float`
Coupling constant

"""

def __init__(self, u, gamma):
super().__init__(value=u, gamma=gamma)
def __init__(self, **kwargs):
if "oid" not in kwargs:
kwargs["value"] = kwargs.pop("u")
super().__init__(**kwargs)

@property
def u(self):
Expand All @@ -547,11 +587,17 @@ class ElectricPotential(_Interpolated):

where :math:`q` is the charge of the particle.

Arguments
----------
field : (M, N, O, 1) array_like of :obj:`float`
Field velocity on a grid of size (M, N, O)
grid_spacing : (3,) array_like of :obj:`float`
Spacing of the grid points.

"""

def __init__(self, field, **kwargs):
super().__init__(field, **kwargs)
def __init__(self, **kwargs):
super().__init__(**kwargs)

_codim = 1
_so_name = "Constraints::ElectricPotential"
23 changes: 22 additions & 1 deletion src/python/espressomd/system.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -75,13 +75,34 @@ if OIF_GLOBAL_FORCES:
cdef bool _system_created = False

cdef class System:
""" The base class for espressomd.system.System().
"""The espresso system class.

.. note:: every attribute has to be declared at the class level.
This means that methods cannot define an attribute by using
``self.new_attr = somevalue`` without declaring it inside this
indentation level, either as method, property or reference.

Attributes
----------
globals : :class:`espressomd.globals.Globals`
actors : :class:`espressomd.actors.Actors`
analysis : :class:`espressomd.analyze.Analysis`
auto_update_accumulators : :class:`espressomd.accumulators.AutoUpdateAccumulators`
bonded_inter : :class:`espressomd.interactions.BondedInteractions`
cell_system : :class:`espressomd.cellsystem.CellSystem`
collision_detection : :class:`espressomd.collision_detection.CollisionDetection`
comfixed : :class:`espressomd.comfixed.ComFixed`
constraints : :class:`espressomd.constraints.Constraints`
cuda_init_handle : :class:`espressomd.cuda_init.CudaInitHandle`
galilei : :class:`espressomd.galilei.GalileiTransform`
integrator : :class:`espressomd.integrate.Integrator`
lbboundaries : :class:`espressomd.lbboundaries.LBBoundaries`
ekboundaries : :class:`espressomd.ekboundaries.EKBoundaries`
minimize_energy : :class:`espressomd.minimize_energy.MinimizeEnergy`
non_bonded_inter : :class:`espressomd.interactions.NonBondedInteractions`
part : :class:`espressomd.particle_data.ParticleList`
thermostat : :class:`espressomd.thermostat.Thermostat`

"""
cdef public:
globals
Expand Down
9 changes: 2 additions & 7 deletions src/script_interface/constraints/fields.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -141,13 +141,8 @@ struct field_params_impl<Interpolated<T, codim>> {
[this_]() { return this_().shape(); }},
{"_field_codim", AutoParameter::read_only,
[]() { return static_cast<int>(codim); }},
{"_field_data", AutoParameter::read_only, [this_]() {
auto &field_data = this_().field_data();
auto data_ptr =
reinterpret_cast<const double *>(field_data.data());
return std::vector<double>(
data_ptr, data_ptr + codim * field_data.num_elements());
}}};
{"_field_data", AutoParameter::read_only,
[this_]() { return this_().field_data_flat(); }}};
}
};

Expand Down
8 changes: 8 additions & 0 deletions src/script_interface/get_value.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,14 @@ struct vector_conversion_visitor : boost::static_visitor<Utils::Vector<T, N>> {
return ret;
}

Utils::Vector<T, N>
operator()(std::vector<T, std::allocator<T>> const &v) const {
if (N != v.size()) {
throw boost::bad_get{};
}
return Utils::Vector<T, N>(v);
}

template <typename U> Utils::Vector<T, N> operator()(U const &) const {
throw boost::bad_get{};
}
Expand Down
Loading