18 for (
Size k = 0; k < f.
size(); ++k) {
40 for (
Size k = 0; k < f.
size(); ++k) {
41 const Size j = neighs[k];
42 const Float u_ji = u[j] - u[i];
43 f[k] = 0.5_f * (1._f + u_ji *
sgn(e[k]) / (
abs(u_ji) + 1._f / (1._f +
abs(u_ji))));
63 for (
Size k = 0; k < f.
size(); ++k) {
64 const Size j = neighs[k];
65 const Float u_ji = u[j] - u[i];
67 const Float A = e[k] / u_ji;
68 const Float B = (A >= 0._f) ? A / m[i] : A / m[j];
74 }
else if (e[k] == 0._f) {
80 f[k] = m[i] / e[k] * ((e[k] + m[i] * u[i] + m[j] * u[j]) / (m[i] + m[j]) - u[i]);
86 template <
typename TPrimary,
typename TSecondary>
96 primary.initialize(storage);
97 secondary.initialize(storage);
106 for (
Size k = 0; k < f.
size(); ++k) {
107 const Size j = neighs[k];
114 f[k] =
lerp(f1, f2, chi);
115 SPH_ASSERT(f[k] >= 0._f && f[k] <= 1._f, f[k]);
128 , threadData(scheduler) {
132 makeAuto<BlendingPartitioner<SmoothlyDiminishingPartitioner, MonotonicDiminishingPartitioner>>();
143 throw InvalidSetup(
"ECS does not support boundary conditions yet");
162 gradList.resize(r.
size());
164 auto evalDerivatives = [&](
const Size i, ThreadData& data) {
167 neighList[i].
clear();
170 for (
auto& n : data.neighs) {
171 const Size j = n.index;
172 const Float hbar = 0.5_f * (r[i][
H] + r[j][
H]);
178 const Vector gr = symmetrizedKernel.grad(r[i], r[j]);
181 neighList[i].
push(j);
182 gradList[i].push(gr);
185 derivatives.
eval(i, neighList[i], gradList[i]);
195 for (ThreadData& data : threadData) {
196 data.energyChange.resize(particleCnt);
197 data.energyChange.fill(0._f);
205 accumulated.
store(storage);
225 auto evalAccelerations = [&](
const Size i, ThreadData& data) {
226 data.accelerations.
resize(neighList[i].size());
229 data.energyChange.resize(neighList[i].size());
230 for (
Size k = 0; k < neighList[i].
size(); ++k) {
231 const Size j = neighList[i][k];
232 const Vector vi12 = v[i] + 0.5_f * dv[i] * dt;
233 const Vector vj12 = v[j] + 0.5_f * dv[j] * dt;
235 data.energyChange[k] = m[i] *
dot(vj12 - vi12, data.accelerations[k]);
238 data.partitions.resize(neighList[i].size());
239 partitioner->
compute(i, neighList[i], data.energyChange, data.partitions);
242 for (
Size k = 0; k < neighList[i].
size(); ++k) {
243 du[i] += data.partitions[k] * data.energyChange[k] / m[i];
250 void EnergyConservingSolver::sanityCheck(
const Storage&
UNUSED(storage))
const {}
INLINE bool isReal(const AntisymmetricTensor &t)
INLINE ArrayView< T > getSingleValueView(T &value)
#define SPH_ASSERT(x,...)
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.
Logging routines of the run.
constexpr INLINE T max(const T &f1, const T &f2)
INLINE T lerp(const T v1, const T v2, const TAmount amount)
constexpr INLINE T sqr(const T &f) noexcept
Return a squared value.
INLINE int sgn(const T val)
INLINE auto abs(const T &f)
#define NAMESPACE_SPH_END
#define MEASURE_SCOPE(name)
@ POSITION
Positions (velocities, accelerations) of particles, always a vector quantity,.
@ ENERGY
Specific internal energy, always a scalar quantity.
@ MASS
Paricles masses, always a scalar quantity.
INLINE void parallelFor(IScheduler &scheduler, const Size from, const Size to, TFunctor &&functor)
Executes a functor concurrently from all available threads.
StaticArray< T0 &, sizeof...(TArgs)+1 > tie(T0 &t0, TArgs &... rest)
Creates a static array from a list of l-value references.
@ RUN_TIME
Current time of the simulation in code units. Does not necessarily have to be 0 when run starts.
@ TIMESTEP_VALUE
Current value of timestep.
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...
void evalAccelerations(const Size idx, ArrayView< const Size > neighs, ArrayView< const Vector > grads, Array< Vector > &dv) const
Storage for accumulating derivatives.
void store(Storage &storage)
Stores accumulated values to corresponding quantities.
INLINE TCounter size() const
void resize(const TCounter newSize)
Resizes the array to new size.
INLINE void push(U &&u)
Adds new element to the end of the array, resizing the array if necessary.
void clear()
Removes all elements from the array, but does NOT release the memory.
INLINE TCounter size() const noexcept
virtual void compute(const Size i, ArrayView< const Size > neighs, ArrayView< const Float > e, ArrayView< Float > f) const override
virtual void initialize(const Storage &storage) override
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.
See Owen 2009: A compatibly differenced total energy conserving form of SPH.
EnergyConservingSolver(IScheduler &scheduler, const RunSettings &settings, const EquationHolder &eqs)
Container holding equation terms.
void setDerivatives(DerivativeHolder &derivatives, const RunSettings &settings) const
Calls EquationTerm::setDerivatives 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.
virtual void compute(const Size UNUSED(i), ArrayView< const Size > UNUSED(neighs), ArrayView< const Float > UNUSED(e), ArrayView< Float > f) const override
virtual void initialize(const Storage &UNUSED(storage)) override
Base class for asymmetric SPH solvers.
IScheduler & scheduler
Scheduler used to parallelize the solver.
Float getMaxSearchRadius(const Storage &storage) const
LutKernel< DIMENSIONS > kernel
Selected SPH kernel.
EquationHolder equations
Holds all equation terms evaluated by the solver.
AutoPtr< ISymmetricFinder > finder
Structure used to search for neighbouring particles.
virtual RawPtr< const IBasicFinder > getFinder(ArrayView< const Vector > r)
Returns a finder, already build using the provided positions.
Interface of objects finding neighbouring particles.
virtual void initialize(const Storage &storage)=0
virtual void compute(const Size i, ArrayView< const Size > neighs, ArrayView< const Float > e, ArrayView< Float > f) const =0
Interface that allows unified implementation of sequential and parallelized versions of algorithms.
Thrown when components of the run are mutually incompatible.
INLINE Float radius() const noexcept
virtual void compute(const Size i, ArrayView< const Size > neighs, ArrayView< const Float > e, ArrayView< Float > f) const override
virtual void initialize(const Storage &storage) override
TValue get(const TEnum idx, std::enable_if_t<!std::is_enum< std::decay_t< TValue >>::value, int >=0) const
Returns a value of given type from the settings.
virtual void compute(const Size i, ArrayView< const Size > neighs, ArrayView< const Float > e, ArrayView< Float > f) const override
virtual void initialize(const Storage &storage) override
Object holding various statistics about current run.
TValue getOr(const StatisticsId idx, const TValue &other) const
Returns value of a statistic, or a given value if the statistic is not stored.
Container storing all quantities used within the simulations.
StaticArray< Array< TValue > &, 3 > getAll(const QuantityId key)
Retrieves quantity buffers from the storage, given its key and value type.
Array< TValue > & getDt(const QuantityId key)
Retrieves a quantity derivative from the storage, given its key and value type.
Size getParticleCnt() const
Returns the number of particles.
Array< TValue > & getValue(const QuantityId key)
Retrieves a quantity values from the storage, given its key and value type.
@ TIMESTEPPING_INITIAL_TIMESTEP