DerivativeHelpers.h

Go to the documentation of this file.

#pragma once
 
#include "quantities/Quantity.h"
#include "sph/equations/Derivative.h"
 
NAMESPACE_SPH_BEGIN
 
enum class DerivativeFlag {
    CORRECTED = 1 << 0,
 
    SUM_ONLY_UNDAMAGED = 1 << 1,
};
 
template <typename TDerived>
class DerivativeTemplate : public ISymmetricDerivative {
private:
    ArrayView<const Size> idxs;
    ArrayView<const Float> reduce;
    ArrayView<const SymmetricTensor> C;
 
    Flags<DerivativeFlag> flags;
 
public:
    explicit DerivativeTemplate(const RunSettings& settings, const Flags<DerivativeFlag> flags = EMPTY_FLAGS)
        : flags(flags) {
        const bool useCorrectionTensor = settings.get<bool>(RunSettingsId::SPH_STRAIN_RATE_CORRECTION_TENSOR);
        if (!useCorrectionTensor) {
            // 'global' override for correction tensor
            this->flags.unset(DerivativeFlag::CORRECTED);
        }
        const bool sumOnlyUndamaged = settings.get<bool>(RunSettingsId::SPH_SUM_ONLY_UNDAMAGED);
        if (!sumOnlyUndamaged) {
            // 'global' override - always sum all particles
            this->flags.unset(DerivativeFlag::SUM_ONLY_UNDAMAGED);
        }
    }
 
    virtual void create(Accumulated& results) override final {
        derived()->additionalCreate(results);
    }
 
    virtual void initialize(const Storage& input, Accumulated& results) override final {
        if (flags.has(DerivativeFlag::CORRECTED)) {
            C = results.getBuffer<SymmetricTensor>(
                QuantityId::STRAIN_RATE_CORRECTION_TENSOR, OrderEnum::ZERO);
        } else {
            C = nullptr;
        }
        if (flags.has(DerivativeFlag::SUM_ONLY_UNDAMAGED) && input.has(QuantityId::STRESS_REDUCING)) {
            idxs = input.getValue<Size>(QuantityId::FLAG);
            reduce = input.getValue<Float>(QuantityId::STRESS_REDUCING);
        } else {
            idxs = nullptr;
            reduce = nullptr;
        }
 
        derived()->additionalInitialize(input, results);
    }
 
    virtual bool equals(const IDerivative& other) const override final {
        if (typeid(*this) != typeid(other)) {
            return false;
        }
        const TDerived* actOther = assert_cast<const TDerived*>(&other);
        return (flags == actOther->flags) && derived()->additionalEquals(*actOther);
    }
 
    virtual void evalNeighs(const Size idx,
        ArrayView<const Size> neighs,
        ArrayView<const Vector> grads) override {
        SPH_ASSERT(neighs.size() == grads.size());
        if (C) {
            sum(idx, neighs, grads, [this](Size i, Size j, const Vector& grad) INL {
                SPH_ASSERT(C[i] != SymmetricTensor::null());
                derived()->template eval<false>(i, j, C[i] * grad);
            });
        } else {
            sum(idx, neighs, grads, [this](Size i, Size j, const Vector& grad) INL {
                derived()->template eval<false>(i, j, grad);
            });
        }
    }
 
    virtual void evalSymmetric(const Size idx,
        ArrayView<const Size> neighs,
        ArrayView<const Vector> grads) override {
        SPH_ASSERT(neighs.size() == grads.size());
        SPH_ASSERT(!flags.has(DerivativeFlag::CORRECTED));
        sum(idx, neighs, grads, [this](Size i, Size j, const Vector& grad) INL {
            derived()->template eval<true>(i, j, grad);
        });
    }
 
private:
    template <typename TFunctor>
    INLINE void sum(const Size i,
        ArrayView<const Size> neighs,
        ArrayView<const Vector> grads,
        const TFunctor& functor) {
        if (reduce) {
            for (Size k = 0; k < neighs.size(); ++k) {
                const Size j = neighs[k];
                if (idxs[i] != idxs[j] || reduce[i] == 0._f || reduce[j] == 0._f) {
                    continue;
                }
                functor(i, j, grads[k]);
            }
        } else {
            for (Size k = 0; k < neighs.size(); ++k) {
                const Size j = neighs[k];
                functor(i, j, grads[k]);
            }
        }
    }
 
    TDerived* derived() {
        return static_cast<TDerived*>(this);
    }
 
    const TDerived* derived() const {
        return static_cast<const TDerived*>(this);
    }
};
 
 
template <typename TDerived>
class AccelerationTemplate : public IAcceleration {
private:
    ArrayView<Vector> dv;
    ArrayView<Float> du;
    ArrayView<const Float> m;
    ArrayView<const Size> idxs;
    ArrayView<const Float> reduce;
    bool sumOnlyUndamaged;
 
public:
    explicit AccelerationTemplate(const RunSettings& settings,
        const Flags<DerivativeFlag> flags = EMPTY_FLAGS) {
        SPH_ASSERT(!flags.has(DerivativeFlag::CORRECTED), "Forces should be never corrected");
 
        // sum only undamaged if specified by the flag and it is specified by the 'global' override
        sumOnlyUndamaged = flags.has(DerivativeFlag::SUM_ONLY_UNDAMAGED) &&
                           settings.get<bool>(RunSettingsId::SPH_SUM_ONLY_UNDAMAGED);
    }
 
    virtual void create(Accumulated& results) override final {
        results.insert<Vector>(QuantityId::POSITION, OrderEnum::SECOND, BufferSource::SHARED);
        results.insert<Float>(QuantityId::ENERGY, OrderEnum::FIRST, BufferSource::SHARED);
 
        derived()->additionalCreate(results);
    }
 
    virtual void initialize(const Storage& input, Accumulated& results) override final {
        dv = results.getBuffer<Vector>(QuantityId::POSITION, OrderEnum::SECOND);
        du = results.getBuffer<Float>(QuantityId::ENERGY, OrderEnum::FIRST);
 
        m = input.getValue<Float>(QuantityId::MASS);
        if (sumOnlyUndamaged && input.has(QuantityId::STRESS_REDUCING)) {
            idxs = input.getValue<Size>(QuantityId::FLAG);
            reduce = input.getValue<Float>(QuantityId::STRESS_REDUCING);
        } else {
            idxs = nullptr;
            reduce = nullptr;
        }
 
        derived()->additionalInitialize(input, results);
    }
 
    virtual bool equals(const IDerivative& other) const override final {
        if (typeid(*this) != typeid(other)) {
            return false;
        }
        const TDerived* actOther = assert_cast<const TDerived*>(&other);
        return (sumOnlyUndamaged == actOther->sumOnlyUndamaged) && derived()->additionalEquals(*actOther);
    }
 
    virtual void evalNeighs(const Size idx,
        ArrayView<const Size> neighs,
        ArrayView<const Vector> grads) override {
        SPH_ASSERT(neighs.size() == grads.size());
        sum(idx, neighs, grads, [this](Size UNUSED(k), Size i, Size j, const Vector& grad) INL {
            Vector f;
            Float de;
            tieToTuple(f, de) = derived()->template eval<false>(i, j, grad);
            dv[i] += m[j] * f;
            du[i] += m[j] * de;
        });
    }
 
    virtual void evalSymmetric(const Size idx,
        ArrayView<const Size> neighs,
        ArrayView<const Vector> grads) override {
        SPH_ASSERT(neighs.size() == grads.size());
        sum(idx, neighs, grads, [this](Size UNUSED(k), Size i, Size j, const Vector& grad) INL {
            Vector f;
            Float de;
            tieToTuple(f, de) = derived()->template eval<true>(i, j, grad);
            dv[i] += m[j] * f;
            dv[j] -= m[i] * f;
            du[i] += m[j] * de;
            du[j] += m[i] * de;
        });
    }
 
    virtual void evalAcceleration(const Size idx,
        ArrayView<const Size> neighs,
        ArrayView<const Vector> grads,
        ArrayView<Vector> dv) override {
        SPH_ASSERT(neighs.size() == grads.size() && neighs.size() == dv.size());
        sum(idx, neighs, grads, [this, &dv](Size k, Size i, Size j, const Vector& grad) INL {
            Vector f;
            tieToTuple(f, IGNORE) = derived()->template eval<false>(i, j, grad);
            dv[k] += m[j] * f;
        });
    }
 
private:
    template <typename TFunctor>
    INLINE void sum(const Size i,
        ArrayView<const Size> neighs,
        ArrayView<const Vector> grads,
        const TFunctor& functor) {
        if (reduce) {
            for (Size k = 0; k < neighs.size(); ++k) {
                const Size j = neighs[k];
                if (idxs[i] != idxs[j] || reduce[i] == 0._f || reduce[j] == 0._f) {
                    continue;
                }
                functor(k, i, j, grads[k]);
            }
        } else {
            for (Size k = 0; k < neighs.size(); ++k) {
                const Size j = neighs[k];
                functor(k, i, j, grads[k]);
            }
        }
    }
 
    INLINE TDerived* derived() {
        return static_cast<TDerived*>(this);
    }
 
    INLINE const TDerived* derived() const {
        return static_cast<const TDerived*>(this);
    }
};
 
class CenterDensityDiscr {
    ArrayView<const Float> rho, m;
 
public:
    void initialize(const Storage& input) {
        tie(rho, m) = input.getValues<Float>(QuantityId::DENSITY, QuantityId::MASS);
    }
 
    template <typename T>
    INLINE T eval(const Size i, const Size j, const T& value) {
        return m[j] / rho[i] * value;
    }
};
 
class NeighbourDensityDiscr {
    ArrayView<const Float> rho, m;
 
public:
    void initialize(const Storage& input) {
        tie(rho, m) = input.getValues<Float>(QuantityId::DENSITY, QuantityId::MASS);
    }
 
    template <typename T>
    INLINE T eval(const Size UNUSED(i), const Size j, const T& value) {
        return m[j] / rho[j] * value;
    }
};
 
 
template <QuantityId Id, typename Discr, typename Traits>
class VelocityTemplate : public DerivativeTemplate<VelocityTemplate<Id, Discr, Traits>> {
private:
    ArrayView<const Vector> v;
 
    Discr discr;
 
    ArrayView<typename Traits::Type> deriv;
 
public:
    using DerivativeTemplate<VelocityTemplate<Id, Discr, Traits>>::DerivativeTemplate;
 
    INLINE void additionalCreate(Accumulated& results) {
        results.insert<typename Traits::Type>(Id, OrderEnum::ZERO, BufferSource::UNIQUE);
    }
 
    INLINE void additionalInitialize(const Storage& input, Accumulated& results) {
        v = input.getDt<Vector>(QuantityId::POSITION);
        discr.initialize(input);
        deriv = results.getBuffer<typename Traits::Type>(Id, OrderEnum::ZERO);
    }
 
    INLINE bool additionalEquals(const VelocityTemplate& UNUSED(other)) const {
        return true;
    }
 
    template <bool Symmetrize>
    INLINE void eval(const Size i, const Size j, const Vector& grad) {
        const typename Traits::Type dv = Traits::eval(v[j] - v[i], grad);
        deriv[i] += discr.eval(i, j, dv);
        if (Symmetrize) {
            deriv[j] += discr.eval(j, i, dv);
        }
    }
};
 
struct DivergenceTraits {
    using Type = Float;
 
    INLINE static Float eval(const Vector& v, const Vector& grad) {
        return dot(v, grad);
    }
};
 
struct RotationTraits {
    using Type = Vector;
 
    INLINE static Vector eval(const Vector& v, const Vector& grad) {
        // nabla x v
        return cross(grad, v);
    }
};
 
struct GradientTraits {
    using Type = SymmetricTensor;
 
    INLINE static SymmetricTensor eval(const Vector& v, const Vector& grad) {
        return symmetricOuter(v, grad);
    }
};
 
 
template <typename Discr>
using VelocityDivergence = VelocityTemplate<QuantityId::VELOCITY_DIVERGENCE, Discr, DivergenceTraits>;
 
template <typename Discr>
using VelocityRotation = VelocityTemplate<QuantityId::VELOCITY_ROTATION, Discr, RotationTraits>;
 
template <typename Discr>
using VelocityGradient = VelocityTemplate<QuantityId::VELOCITY_GRADIENT, Discr, GradientTraits>;
 
 
template <template <class> class TVelocityDerivative>
AutoPtr<IDerivative> makeDerivative(const RunSettings& settings, Flags<DerivativeFlag> flags = EMPTY_FLAGS) {
    const DiscretizationEnum formulation =
        settings.get<DiscretizationEnum>(RunSettingsId::SPH_DISCRETIZATION);
    switch (formulation) {
    case DiscretizationEnum::STANDARD:
        return makeAuto<TVelocityDerivative<CenterDensityDiscr>>(settings, flags);
    case DiscretizationEnum::BENZ_ASPHAUG:
        return makeAuto<TVelocityDerivative<NeighbourDensityDiscr>>(settings, flags);
    default:
        NOT_IMPLEMENTED;
    }
}
 
NAMESPACE_SPH_END