XSph.h

Go to the documentation of this file.

#pragma once
 
 
#include "sph/equations/DerivativeHelpers.h"
#include "sph/equations/EquationTerm.h"
#include "sph/kernel/Kernel.h"
#include "system/Factory.h"
 
NAMESPACE_SPH_BEGIN
 
 
class XSph : public IEquationTerm {
private:
    class Derivative : public DerivativeTemplate<Derivative> {
    private:
        SymmetrizeSmoothingLengths<LutKernel<DIMENSIONS>> kernel;
 
        ArrayView<const Vector> r, v;
        ArrayView<const Float> rho, m;
        ArrayView<Vector> dr;
        Float epsilon;
 
    public:
        explicit Derivative(const RunSettings& settings)
            : DerivativeTemplate<Derivative>(settings)
            , kernel(Factory::getKernel<DIMENSIONS>(settings)) {
            epsilon = settings.get<Float>(RunSettingsId::SPH_XSPH_EPSILON);
        }
 
        INLINE void additionalCreate(Accumulated& results) {
            results.insert<Vector>(QuantityId::XSPH_VELOCITIES, OrderEnum::ZERO, BufferSource::UNIQUE);
        }
 
        INLINE void additionalInitialize(const Storage& input, Accumulated& results) {
            dr = results.getBuffer<Vector>(QuantityId::XSPH_VELOCITIES, OrderEnum::ZERO);
            tie(rho, m) = input.getValues<Float>(QuantityId::DENSITY, QuantityId::MASS);
            ArrayView<const Vector> dummy;
            tie(r, v, dummy) = input.getAll<Vector>(QuantityId::POSITION);
        }
 
        INLINE bool additionalEquals(const Derivative& other) const {
            return epsilon == other.epsilon;
        }
 
        template <bool Symmetric>
        INLINE void eval(const Size i, const Size j, const Vector& UNUSED(grad)) {
            // this depends on v[i]-v[j], so it is zero for i==j
            Vector f = epsilon * (v[j] - v[i]) / (0.5_f * (rho[i] + rho[j])) * kernel.value(r[i], r[j]);
            f[H] = 0._f;
            dr[i] += m[j] * f;
            if (Symmetric) {
                dr[j] -= m[i] * f;
            }
        }
    };
 
public:
    virtual void setDerivatives(DerivativeHolder& derivatives, const RunSettings& settings) override {
        derivatives.require(makeAuto<Derivative>(settings));
    }
 
    virtual void initialize(IScheduler& UNUSED(scheduler), Storage& storage, const Float UNUSED(t)) override {
        // fix previously modified velocities before computing derivatives
        ArrayView<Vector> v = storage.getDt<Vector>(QuantityId::POSITION);
        ArrayView<Vector> dr = storage.getValue<Vector>(QuantityId::XSPH_VELOCITIES);
        for (Size i = 0; i < v.size(); ++i) {
            v[i] -= dr[i];
        }
    }
 
    virtual void finalize(IScheduler& UNUSED(scheduler), Storage& storage, const Float UNUSED(t)) override {
        ArrayView<Vector> v = storage.getDt<Vector>(QuantityId::POSITION);
        ArrayView<Vector> dr = storage.getValue<Vector>(QuantityId::XSPH_VELOCITIES);
        for (Size i = 0; i < v.size(); ++i) {
            v[i] += dr[i];
        }
    }
 
    virtual void create(Storage& storage, IMaterial& UNUSED(material)) const override {
        storage.insert<Vector>(QuantityId::XSPH_VELOCITIES, OrderEnum::ZERO, Vector(0._f));
    }
};
 
NAMESPACE_SPH_END