44 template <
bool Symmetrize>
48 dv[i] +=
Y[i] *
Y[j] / m[i] * (1._f / y[i] + 1._f / y[j]) * grad;
51 du[i] +=
Y[i] *
Y[j] / m[i] / y[i] *
dot(v[i] - v[j], grad);
59 derivatives.
require(makeAuto<DensityIndependentPressureGradient>(settings));
65 for (
Size i = 0; i < p.
size(); ++i) {
67 p[i] =
max(p[i], 100._f);
77 throw InvalidSetup(
"DISPH requires EosMaterial or derived");
86 : scheduler(scheduler)
87 , threadData(scheduler) {
89 kernel = Factory::getKernel<3>(settings);
93 equations += makeTerm<DensityIndependentPressureForce>();
94 equations += makeTerm<SolidStressForce>(settings);
95 equations += makeTerm<StandardAV>();
104 finder->
build(scheduler, r);
112 for (
Size i = 0; i < r.
size(); ++i) {
113 rho[i] = m[i] * y[i] /
Y[i];
114 SPH_ASSERT(rho[i] > 1._f && rho[i] < 1.e4_f, rho[i]);
120 material->initialize(scheduler, storage, material.
sequence());
129 for (
Size i = 0; i < r.
size(); ++i) {
130 Y[i] = m[i] * p[i] / rho[i];
137 auto pressureFunc = [
this,
Y, r, p, &y,
radius](
const Size i, ThreadData& data) {
144 y[i] +=
Y[j] * kernel.
value(r[i], r[j]);
154 auto equationFunc = [
this,
Y, r,
radius](
const Size i, ThreadData& data) {
163 const Float hbar = 0.5_f * (r[i][
H] + r[j][
H]);
172 data.grads.emplaceBack(gr);
173 data.idxs.emplaceBack(j);
176 derivatives.
eval(i, data.idxs, data.grads);
181 equations.
finalize(scheduler, storage, t);
188 for (
Size i = 0; i < r.
size(); ++i) {
189 const Float gamma = rho[i] / p[i] *
sqr(cs[i]);
193 dY[i] = 0._f * ( gamma - 1._f) * m[i] * du[i];
199 material->finalize(scheduler, storage, material.
sequence());
212 equations.
create(storage, material);
@ SHARED
Multiple derivatives may accumulate into the buffer.
INLINE bool isReal(const AntisymmetricTensor &t)
#define SPH_ASSERT(x,...)
Density-independent formulation of SPH.
uint32_t Size
Integral type used to index arrays (by default).
double Float
Precision used withing the code. Use Float instead of float or double where precision is important.
SPH-specific implementation of particle materials.
constexpr INLINE T max(const T &f1, const T &f2)
constexpr INLINE T sqr(const T &f) noexcept
Return a squared value.
INLINE auto abs(const T &f)
#define INLINE
Macros for conditional compilation based on selected compiler.
#define NAMESPACE_SPH_END
@ PRESSURE
Pressure, affected by yielding and fragmentation model, always a scalar quantity.
@ POSITION
Positions (velocities, accelerations) of particles, always a vector quantity,.
@ ENERGY
Specific internal energy, always a scalar quantity.
@ DENSITY
Density, always a scalar quantity.
@ MASS
Paricles masses, always a scalar quantity.
@ SOUND_SPEED
Sound speed, always a scalar quantity.
@ GENERALIZED_ENERGY
The "Y" quantity defined by , used to compute equation of motion and energy in DISPH.
@ SECOND
Quantity with 1st and 2nd derivative.
@ FIRST
Quantity with 1st derivative.
@ ZERO
Quantity without derivatives, or "zero order" of quantity.
INLINE void parallelFor(IScheduler &scheduler, const Size from, const Size to, TFunctor &&functor)
Executes a functor concurrently from all available threads.
Standard SPH artificial viscosity.
StaticArray< T0 &, sizeof...(TArgs)+1 > tie(T0 &t0, TArgs &... rest)
Creates a static array from a list of l-value references.
Statistics gathered and periodically displayed during the run.
@ RUN_TIME
Current time of the simulation in code units. Does not necessarily have to be 0 when run starts.
INLINE Float getSqrLength(const Vector &v)
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...
Storage for accumulating derivatives.
Array< TValue > & getBuffer(const QuantityId id, const OrderEnum order)
Returns the buffer of given quantity and given order.
void store(Storage &storage)
Stores accumulated values to corresponding quantities.
void insert(const QuantityId id, const OrderEnum order, const BufferSource source)
Creates a new storage with given ID.
Object providing safe access to continuous memory of data.
INLINE TCounter size() const
virtual void finalize(IScheduler &UNUSED(scheduler), Storage &UNUSED(storage), const Float UNUSED(t)) override
virtual void initialize(IScheduler &UNUSED(scheduler), Storage &storage, const Float UNUSED(t)) override
virtual void create(Storage &storage, IMaterial &material) const override
Creates all quantities needed by the term using given material.
virtual void setDerivatives(DerivativeHolder &derivatives, const RunSettings &settings) override
Sets derivatives required by this term.
INLINE void additionalCreate(Accumulated &results)
INLINE void additionalInitialize(const Storage &input, Accumulated &results)
DensityIndependentPressureGradient(const RunSettings &settings)
INLINE bool additionalEquals(const DensityIndependentPressureGradient &UNUSED(other)) const
INLINE void eval(const Size i, const Size j, const Vector &grad)
DensityIndependentSolver(IScheduler &scheduler, const RunSettings &settings)
virtual void integrate(Storage &storage, Statistics &stats) override
Computes derivatives of all time-dependent quantities.
~DensityIndependentSolver()
virtual void create(Storage &storage, IMaterial &material) const override
Initializes all quantities needed by the solver in the storage.
Container holding derivatives and the storage they accumulate to.
virtual void require(AutoPtr< IDerivative > &&derivative)
Adds derivative if not already present.
void eval(const Size idx, ArrayView< const Size > neighs, ArrayView< const Vector > grads)
Evaluates all held derivatives for given particle.
INLINE Accumulated & getAccumulated()
virtual void initialize(const Storage &input)
Initialize derivatives before loop.
Helper template for derivatives that define both the symmetrized and asymmetric variant.
Material holding equation of state.
void setDerivatives(DerivativeHolder &derivatives, const RunSettings &settings) const
Calls EquationTerm::setDerivatives for all stored equation terms.
void create(Storage &storage, IMaterial &material) const
Calls EquationTerm::create for all stored equation terms.
void initialize(IScheduler &scheduler, Storage &storage, const Float t)
Calls EquationTerm::initialize for all stored equation terms.
void finalize(IScheduler &scheduler, Storage &storage, const Float t)
Calls EquationTerm::finalize for all stored equation terms.
void build(IScheduler &scheduler, ArrayView< const Vector > points)
Constructs the struct with an array of vectors.
virtual Size findAll(const Size index, const Float radius, Array< NeighbourRecord > &neighbours) const =0
Finds all neighbours within given radius from the point given by index.
Represents a term or terms appearing in evolutionary equations.
Material settings and functions specific for one material.
INLINE TValue getParam(const BodySettingsId paramIdx) const
Returns a parameter associated with given particle.
void setRange(const QuantityId id, const Interval &range, const Float minimal)
Sets the timestepping parameters of given quantity.
Interface that allows unified implementation of sequential and parallelized versions of algorithms.
Object representing a 1D interval of real numbers.
Thrown when components of the run are mutually incompatible.
Non-owning wrapper of a material and particles with this material.
INLINE IndexSequence sequence()
Returns iterable index sequence.
Object holding various statistics about current run.
TValue get(const StatisticsId idx) const
Returns value of a statistic.
Container storing all quantities used within the simulations.
Size getMaterialCnt() const
Return the number of materials in the storage.
Quantity & insert(const QuantityId key, const OrderEnum order, const TValue &defaultValue)
Creates a quantity in the storage, given its key, value type and order.
Array< TValue > & getDt(const QuantityId key)
Retrieves a quantity derivative from the storage, given its key and value type.
auto getValues(const QuantityId first, const QuantityId second, const TArgs... others)
Retrieves an array of quantities from the key.
MaterialView getMaterial(const Size matIdx) const
Returns an object containing a reference to given material.
Array< TValue > & getValue(const QuantityId key)
Retrieves a quantity values from the storage, given its key and value type.
INLINE Float value(const Vector &r1, const Vector &r2) const
INLINE Float radius() const
INLINE Vector grad(const Vector &r1, const Vector &r2) const
Creating code components based on values from settings.
@ ENERGY_MIN
Estimated minimal value of energy used to determine timestepping error.
@ ENERGY
Initial specific internal energy.
@ DENSITY
Density at zero pressure.
@ ENERGY_RANGE
Allowed range of specific internal energy.
AutoPtr< ISymmetricFinder > getFinder(const RunSettings &settings)
Holds information about a neighbour particles.
Size index
Index of particle in the storage.
Float distanceSqr
Squared distance of the particle from the queried particle / position.