1 #include "sph/solvers/SummationSolver.h"
2 #include "objects/finders/NeighborFinder.h"
3 #include "quantities/IMaterial.h"
4 #include "sph/boundary/Boundary.h"
5 #include "sph/equations/av/Standard.h"
6 #include "sph/kernel/Kernel.h"
7 #include "system/Factory.h"
8 #include "system/Statistics.h"
9 #include "thread/AtomicFloat.h"
10
11 NAMESPACE_SPH_BEGIN
12
getEquations(const RunSettings & settings)13 static EquationHolder getEquations(const RunSettings& settings) {
14 Flags<ForceEnum> forces = settings.getFlags<ForceEnum>(RunSettingsId::SPH_SOLVER_FORCES);
15 EquationHolder equations;
16 if (forces.has(ForceEnum::PRESSURE)) {
17 equations += makeTerm<PressureForce>();
18 }
19 if (forces.has(ForceEnum::SOLID_STRESS)) {
20 equations += makeTerm<SolidStressForce>(settings);
21 }
22 if (auto av = Factory::getArtificialViscosity(settings)) {
23 equations += EquationHolder(std::move(av));
24 }
25
26 // we evolve density and smoothing length ourselves (outside the equation framework),
27 // so make sure it does not change outside the solver
28 equations += makeTerm<ConstSmoothingLength>();
29
30 SPH_ASSERT(
31 !forces.has(ForceEnum::SELF_GRAVITY), "Summation solver cannot be currently used with gravity");
32
33 return equations;
34 }
35
36 template <Size Dim>
SummationSolver(IScheduler & scheduler,const RunSettings & settings,const EquationHolder & additionalEquations,AutoPtr<IBoundaryCondition> && bc)37 SummationSolver<Dim>::SummationSolver(IScheduler& scheduler,
38 const RunSettings& settings,
39 const EquationHolder& additionalEquations,
40 AutoPtr<IBoundaryCondition>&& bc)
41 : SymmetricSolver<Dim>(scheduler, settings, getEquations(settings) + additionalEquations, std::move(bc)) {
42
43 targetDensityDifference = settings.get<Float>(RunSettingsId::SPH_SUMMATION_DENSITY_DELTA);
44 densityKernel = Factory::getKernel<Dim>(settings);
45 Flags<SmoothingLengthEnum> flags =
46 settings.getFlags<SmoothingLengthEnum>(RunSettingsId::SPH_ADAPTIVE_SMOOTHING_LENGTH);
47 adaptiveH = (flags != EMPTY_FLAGS);
48 maxIteration = adaptiveH ? settings.get<int>(RunSettingsId::SPH_SUMMATION_MAX_ITERATIONS) : 1;
49 }
50
51 template <Size Dim>
SummationSolver(IScheduler & scheduler,const RunSettings & settings,const EquationHolder & additionalEquations)52 SummationSolver<Dim>::SummationSolver(IScheduler& scheduler,
53 const RunSettings& settings,
54 const EquationHolder& additionalEquations)
55 : SummationSolver(scheduler, settings, additionalEquations, Factory::getBoundaryConditions(settings)) {}
56
57 template <Size Dim>
create(Storage & storage,IMaterial & material) const58 void SummationSolver<Dim>::create(Storage& storage, IMaterial& material) const {
59 const Float rho0 = material.getParam<Float>(BodySettingsId::DENSITY);
60 storage.insert<Float>(QuantityId::DENSITY, OrderEnum::ZERO, rho0);
61 material.setRange(QuantityId::DENSITY, BodySettingsId::DENSITY_RANGE, BodySettingsId::DENSITY_MIN);
62 storage.insert<Size>(QuantityId::NEIGHBOR_CNT, OrderEnum::ZERO, 0);
63 this->equations.create(storage, material);
64 }
65
66 template <Size Dim>
beforeLoop(Storage & storage,Statistics & stats)67 void SummationSolver<Dim>::beforeLoop(Storage& storage, Statistics& stats) {
68 SymmetricSolver<Dim>::beforeLoop(storage, stats);
69 ArrayView<Vector> r = storage.getValue<Vector>(QuantityId::POSITION);
70 ArrayView<Float> m = storage.getValue<Float>(QuantityId::MASS);
71
72 // initialize density estimate
73 rho.resize(r.size());
74 rho.fill(EPS);
75
76 // initialize smoothing length to current values
77 h.resize(r.size());
78 for (Size i = 0; i < r.size(); ++i) {
79 h[i] = r[i][H];
80 SPH_ASSERT(h[i] > 0._f);
81 }
82
83 Float eta = 0._f;
84 for (Size matId = 0; matId < storage.getMaterialCnt(); ++matId) {
85 // all eta's should be the same, but let's use max to be sure
86 eta = max(eta, storage.getMaterial(matId)->getParam<Float>(BodySettingsId::SMOOTHING_LENGTH_ETA));
87 }
88
89 Atomic<Float> totalDiff = 0._f; /// \todo use thread local for summing?
90 auto functor = [this, r, m, eta, &totalDiff](const Size i, ThreadData& data) {
91 /// \todo do we have to recompute neighbors in every iteration?
92 // find all neighbors
93 this->finder->findAll(i, h[i] * densityKernel.radius(), data.neighs);
94 SPH_ASSERT(data.neighs.size() > 0, data.neighs.size());
95 // find density and smoothing length by self-consistent solution.
96 const Float rho0 = rho[i];
97 rho[i] = 0._f;
98 for (auto& n : data.neighs) {
99 const Size j = n.index;
100 /// \todo can this be generally different kernel than the one used for derivatives?
101 rho[i] += m[j] * densityKernel.value(r[i] - r[j], h[i]);
102 }
103 SPH_ASSERT(rho[i] > 0._f, rho[i]);
104 h[i] = eta * root<Dim>(m[i] / rho[i]);
105 SPH_ASSERT(h[i] > 0._f);
106 totalDiff += abs(rho[i] - rho0) / (rho[i] + rho0);
107 };
108
109 this->finder->build(this->scheduler, r);
110 Size iterationIdx = 0;
111 for (; iterationIdx < maxIteration; ++iterationIdx) {
112 parallelFor(this->scheduler, this->threadData, 0, r.size(), functor);
113 const Float diff = totalDiff / r.size();
114 if (diff < targetDensityDifference) {
115 break;
116 }
117 }
118 stats.set(StatisticsId::SOLVER_SUMMATION_ITERATIONS, int(iterationIdx));
119 // save computed values
120 std::swap(storage.getValue<Float>(QuantityId::DENSITY), rho);
121 if (adaptiveH) {
122 for (Size i = 0; i < r.size(); ++i) {
123 r[i][H] = h[i];
124 SPH_ASSERT(r[i][H] > 0._f);
125 }
126 }
127 }
128
129 template <Size Dim>
sanityCheck(const Storage & UNUSED (storage)) const130 void SummationSolver<Dim>::sanityCheck(const Storage& UNUSED(storage)) const {
131 // we handle smoothing lengths ourselves, bypass the check of equations
132 }
133
134 template class SummationSolver<1>;
135 template class SummationSolver<2>;
136 template class SummationSolver<3>;
137
138 NAMESPACE_SPH_END
139