SPH
Creating custom forces, heating or cooling terms, etc.

Code allows for adding additional terms to the existing solver. Any equation terms shall be derived from IEquationTerm and implement the virtual functions.

Harmonic oscillator

EquationTerm can be used to add external potential of harmonic oscillator,

class SpringForce : public IEquationTerm {
private:
Float k;
public:
SpringForce(const Float springConstant)
: k(springConstant) {}
virtual void setDerivatives(DerivativeHolder&, const RunSettings&) override {
// No derivatives necessary, see below.
}
virtual void initialize(Storage&) override {
// Initialize the equation term before evaluating the derivatives. Could be used to enforce some
// dependence between quantities, clamp their values, etc. In this case, we don't have to do anything.
}
virtual void finalize(Storage& storage) override {
// Get particle positions, their masses and accelerations from the storage.
ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASSES);
// Add acceleration due to spring
for (Size i = 0; i < r.size(); ++i) {
dv[i] += -k / m[i] * r[i];
}
}
virtual void create(Storage&, Abstract::Material&) const override {
// No additional quantities are needed by this force.
}
};
uint32_t Size
Integral type used to index arrays (by default).
Definition: Globals.h:16
double Float
Precision used withing the code. Use Float instead of float or double where precision is important.
Definition: Globals.h:13
StaticArray< Vector, 8 > POSITIONS
Definition: Presets.cpp:341
Object providing safe access to continuous memory of data.
Definition: ArrayView.h:17
INLINE TCounter size() const
Definition: ArrayView.h:101
Container holding derivatives and the storage they accumulate to.
Definition: Derivative.h:173
Represents a term or terms appearing in evolutionary equations.
Definition: EquationTerm.h:22
virtual void finalize(IScheduler &scheduler, Storage &storage, const Float t)=0
Computes all the derivatives and/or quantity values based on accumulated derivatives.
virtual void setDerivatives(DerivativeHolder &derivatives, const RunSettings &settings)=0
Sets derivatives required by this term.
virtual void initialize(IScheduler &scheduler, Storage &storage, const Float t)=0
Initialize all the derivatives and/or quantity values before derivatives are computed.
virtual void create(Storage &storage, IMaterial &material) const =0
Creates all quantities needed by the term using given material.
Container storing all quantities used within the simulations.
Definition: Storage.h:230
Array< TValue > & getD2t(const QuantityId key)
Retrieves a quantity second derivative from the storage, given its key and value type.
Definition: Storage.cpp:243
Array< TValue > & getValue(const QuantityId key)
Retrieves a quantity values from the storage, given its key and value type.
Definition: Storage.cpp:191

Harmonic force only depends on the position of particle, independent for other particles. In other words, doesn't depend on spatial derivatives of any quantity. PressureForce, on the other hand, does depend. Implementations of IEquationTerm never directly compute interactions between particles. This is handled by different object, IDerivative.

Attention
Equation terms should not modify quantity values directly, nor should they assign values of derivatives. It is only recommended to modify highest-order derivatives by addition or subtraction, as in the example. This ensures the resulting derivatives and subsequently quantity values are independent on order of evaluation of equation terms. Direct modification of quantity can be only done if the quantity does not appear in other equation terms.

Heat diffusion equation

// We can create additional quantities beside those listed in QuantityId enum, preferably with negative
// indices to avoid conflicts between quantities. Any quantity index have to be cast to QuantityId.
enum class AdditionalId {
};
// Object computing laplacian of internal energy by SPH dicretization.
class EnergyLaplacian : public IDerivative {
private:
// View of the buffer where we save computed values
// Input quantities
public:
virtual void create(Accumulated& results) override {
// Create a new buffer for storing the energy laplacian. Note that this is called for every thread;
// every thread computes its own derivatives and these thread-local derivatives are summed up after
// every particle pair has been evaluated.
results.insert<Float>(QuantityId(AdditionalId::ENERGY_LAPLACIAN));
}
virtual void initialize(const Storage& input, Accumulated& results) override {
// Fetch the array views of input quantities and view of the output buffer
tie(u, m, rho) = input.getValues<Float>(QuantityId::ENERGY, QuantityId::MASSES, QuantityId::DENSITY);
flow = results.getValue<Float>(QuantityId(AdditionalId::ENERGY_LAPLACIAN));
}
virtual void evalSymmetric(const Size i, ArrayView<const Size> neighs, ArrayView<const Vector> grads) override {
// Loop over particle neighbours, computes a part of the sum for both of interacting particles. Derivatives
// for both particles have to be computed. ArrayView neighs contains indices of neighbours, grads contains
// precomputed gradients of the SPH kernel.
for (Size k = 0; k < neighs.size(); ++k) {
const Size j = neighs[k];
// SPH laplacian approximation, see Price (2010)
const Float laplacian = 2._f * dot(r[j] - r[i], grads[k]) / getSqrLength(r[j] - r[i]);
const Float f = (u[j] - u[i]) * laplacian;
flow[i] += m[j] / rho[j] * f;
flow[j] -= m[i] / rho[i] * f;
}
}
virtual void evalAsymmetric(const Size i, ArrayView<const Size> neighs, ArrayView<const Vector> grads) override {
// Asymmetric variant of the function, summing up only derivatives for i-th particle. Note that this code
// duplication can be avoided by deriving from DerivativeTemplate rather than directly from
// IDerivative.
for (Size k = 0; k < neighs.size(); ++k) {
const Size j = neighs[k];
const Float laplacian = 2._f * dot(r[j] - r[i], grads[k]) / getSqrLength(r[j] - r[i]);
const Float f = (u[j] - u[i]) * laplacian;
flow[i] += m[j] / rho[j] * f;
}
}
};
// Right-hand side of the heat diffusion equation, uses the above defined derivative to compute
// the final energy derivatives.
private:
Float alpha; // diffusivity
public:
: alpha(alpha) {}
virtual void setDerivatives(DerivativeHolder& derivatives, const RunSettings& settings) override {
// Add laplacian of energy to the list of derivatives we wish to evaluate
derivatives.require<EnergyLaplacian>(settings);
}
virtual void initialize(Storage&) override {
// Nothing to be done here.
}
virtual void finalize(Storage& storage) override {
// When finalize is called, the solver already computed all derivatives, summed up
// the thread-local buffers and stored the results into the particle storage. Here we can
// simply retrieve the computed laplacian and write the final equation for energy derivative.
ArrayView<Float> delta = storage.getValue<Float>(QuantityId(AdditionalId::ENERGY_LAPLACIAN));
for (Size i = 0; i < du.size(); ++i) {
// Heat diffusion equation
du[i] += -alpha * delta[i];
}
}
virtual void create(Storage& storage, Abstract::Material&) const override {
// Create a quantity in the storage for the energy laplacian, to make sure the computed derivatives
// can be stored somewhere. This quantity can be later saved to output file for debugging purposes.
// Alternatively, we could directly increment the energy derivatives by the computed laplacian,
// to avoid allocating unnecessary buffers.
storage.insert<Float>(QuantityId(AdditionalId::ENERGY_LAPLACIAN), OrderEnum::ZERO, 0._f);
}
};
INLINE T laplacian(const T &value, const Vector &grad, const Vector &dr)
SPH approximation of laplacian, computed from a kernel gradient.
Definition: Kernel.h:519
QuantityId
Unique IDs of basic quantities of SPH particles.
Definition: QuantityIds.h:19
@ ENERGY_LAPLACIAN
Laplacian of internal energy, used in heat diffusion equation.
@ ENERGY
Specific internal energy, always a scalar quantity.
@ DENSITY
Density, always a scalar quantity.
@ ZERO
Quantity without derivatives, or "zero order" of quantity.
StaticArray< T0 &, sizeof...(TArgs)+1 > tie(T0 &t0, TArgs &... rest)
Creates a static array from a list of l-value references.
Definition: StaticArray.h:281
INLINE Float getSqrLength(const Vector &v)
Definition: Vector.h:574
INLINE float dot(const BasicVector< float > &v1, const BasicVector< float > &v2)
Make sure the vector is trivially constructible and destructible, needed for fast initialization of a...
Definition: Vector.h:548
Storage for accumulating derivatives.
Definition: Accumulated.h:30
void insert(const QuantityId id, const OrderEnum order, const BufferSource source)
Creates a new storage with given ID.
virtual void require(AutoPtr< IDerivative > &&derivative)
Adds derivative if not already present.
Definition: Derivative.cpp:77
virtual void create(Accumulated &results) override final
virtual void evalSymmetric(const Size idx, ArrayView< const Size > neighs, ArrayView< const Vector > grads) override
virtual void initialize(const Storage &input, Accumulated &results) override final
virtual void initialize(IScheduler &UNUSED(scheduler), Storage &UNUSED(storage), const Float UNUSED(t)) override
Definition: Heat.h:57
virtual void finalize(IScheduler &UNUSED(scheduler), Storage &storage, const Float UNUSED(t)) override
Definition: Heat.h:61
virtual void create(Storage &storage, IMaterial &material) const override
Creates all quantities needed by the term using given material.
Definition: Heat.h:73
virtual void setDerivatives(DerivativeHolder &derivatives, const RunSettings &settings) override
Sets derivatives required by this term.
Definition: Heat.h:52
Derivative accumulated by summing up neighbouring particles.
Definition: Derivative.h:43
Quantity & insert(const QuantityId key, const OrderEnum order, const TValue &defaultValue)
Creates a quantity in the storage, given its key, value type and order.
Definition: Storage.cpp:270
Array< TValue > & getDt(const QuantityId key)
Retrieves a quantity derivative from the storage, given its key and value type.
Definition: Storage.cpp:217
auto getValues(const QuantityId first, const QuantityId second, const TArgs... others)
Retrieves an array of quantities from the key.
Definition: Storage.h:359