ViennaLS
Loading...
Searching...
No Matches
lsStencilLocalLaxFriedrichsScalar.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <hrleSparseBoxIterator.hpp>
4
5#include <lsDomain.hpp>
6#include <lsExpand.hpp>
8#include <lsVelocityField.hpp>
9
10#include <vcVectorType.hpp>
11
12namespace lsInternal {
13
14using namespace viennacore;
15
22template <class T, int D, int order,
23 DifferentiationSchemeEnum finiteDifferenceScheme =
26 using LevelSetType = SmartPointer<viennals::Domain<T, D>>;
27 using LevelSetsType = std::vector<LevelSetType>;
28
29 LevelSetType levelSet;
30 SmartPointer<viennals::VelocityField<T>> velocities;
31 viennahrle::SparseBoxIterator<viennahrle::Domain<T, D>,
32 static_cast<int>(finiteDifferenceScheme) + 1 +
33 order>
34 neighborIterator;
35 const double alphaFactor;
36 const double normalEpsilon =
37 std::cbrt(std::numeric_limits<double>::epsilon());
38
39 // Final dissipation coefficients that are used by the time integrator. If
40 // D==2 last entries are 0.
41 Vec3D<T> finalAlphas;
42 static constexpr unsigned numStencilPoints = hrleUtil::pow(2 * order + 1, D);
43
44 static T pow2(const T &value) { return value * value; }
45
46 Vec3D<T> calculateNormal(const viennahrle::Index<D> &offset) {
47 Vec3D<T> normal = {0.0, 0.0, 0.0};
48 constexpr int startIndex = -1;
49 T modulus = 0.;
50
51 for (unsigned i = 0; i < D; ++i) {
52 viennahrle::Index<D> index(offset);
53 std::vector<T> values;
54 for (unsigned j = 0; j < 3; ++j) {
55 index[i] = startIndex + j;
56 values.push_back(neighborIterator.getNeighbor(index).getValue());
57 }
59 values.data(), levelSet->getGrid().getGridDelta());
60 modulus += normal[i] * normal[i];
61 }
62 modulus = std::sqrt(modulus);
63 for (unsigned i = 0; i < D; ++i) {
64 normal[i] /= modulus;
65 }
66 return normal;
67 }
68
69 VectorType<T, D> calculateGradient(const viennahrle::Index<D> &offset) {
70 VectorType<T, D> gradient;
71
72 constexpr unsigned numValues =
74 const int startIndex = -std::floor(numValues / 2);
75
76 for (unsigned i = 0; i < D; ++i) {
77 viennahrle::Index<D> index(offset);
78 std::vector<T> values;
79 for (unsigned j = 0; j < numValues; ++j) {
80 index[i] = startIndex + j;
81 values.push_back(neighborIterator.getNeighbor(index).getValue());
82 }
83
84 gradient[i] =
86 values.data(), levelSet->getGrid().getGridDelta());
87 }
88
89 return gradient;
90 }
91
92 VectorType<T, D> calculateGradientDiff() {
93 VectorType<T, D> gradient;
94
95 const unsigned numValues =
97 const int startIndex = -std::floor(numValues / 2);
98
99 for (unsigned i = 0; i < D; ++i) {
100 viennahrle::Index<D> index(0);
101 std::vector<T> values;
102 for (unsigned j = 0; j < numValues; ++j) {
103 index[i] = startIndex + j;
104 values.push_back(neighborIterator.getNeighbor(index).getValue());
105 }
106
107 gradient[i] =
109 &(values[0]), levelSet->getGrid().getGridDelta());
110 }
111
112 return gradient;
113 }
114
115public:
116 static void prepareLS(LevelSetType passedlsDomain) {
117 // Expansion of sparse field must depend on spatial derivative order
118 // AND slf stencil order! --> currently assume scheme = 3rd order always
119 viennals::Expand<T, D>(passedlsDomain, 2 * (order + 1) + 4).apply();
120 }
121
122 StencilLocalLaxFriedrichsScalar(LevelSetType passedlsDomain,
123 SmartPointer<viennals::VelocityField<T>> vel,
124 double a = 1.0)
125 : levelSet(passedlsDomain), velocities(vel),
126 neighborIterator(levelSet->getDomain()), alphaFactor(a) {
127 for (int i = 0; i < 3; ++i) {
128 finalAlphas[i] = 0;
129 }
130 }
131
132 std::pair<T, T> operator()(const viennahrle::Index<D> &indices,
133 int material) {
134 auto &grid = levelSet->getGrid();
135 double gridDelta = grid.getGridDelta();
136
137 Vec3D<T> coordinate{0., 0., 0.};
138 for (unsigned i = 0; i < D; ++i) {
139 coordinate[i] = indices[i] * gridDelta;
140 }
141
142 // move neighborIterator to current position
143 neighborIterator.goToIndicesSequential(indices);
144
145 // convert coordinate to std array for interface
146 Vec3D<T> coordArray{coordinate[0], coordinate[1], coordinate[2]};
147
148 // if there is a vector velocity, we need to project it onto a scalar
149 // velocity first using its normal vector
150 // /*if (vectorVelocity != Vec3D<T>({}))*/ {
151 Vec3D<T> normalVector;
152 T denominator = 0; // normal modulus
153 for (unsigned i = 0; i < D; i++) {
154 viennahrle::Index<D> neighborIndex(0);
155 neighborIndex[i] = 1;
156 // normal vector calculation
157 T pos = neighborIterator.getNeighbor(neighborIndex).getValue() -
158 neighborIterator.getCenter().getValue();
159 T neg = neighborIterator.getCenter().getValue() -
160 neighborIterator.getNeighbor(-neighborIndex).getValue();
161 normalVector[i] = (pos + neg) * 0.5;
162 // normalise normal vector
163 denominator += normalVector[i] * normalVector[i];
164 }
165 denominator = std::sqrt(denominator);
166
167 for (unsigned i = 0; i < D; ++i) {
168 normalVector[i] /= denominator;
169 }
170
171 double scalarVelocity = velocities->getScalarVelocity(
172 coordArray, material, normalVector,
173 neighborIterator.getCenter().getPointId());
174 auto vectorVelocity = velocities->getVectorVelocity(
175 coordArray, material, normalVector,
176 neighborIterator.getCenter().getPointId());
177
178 // now calculate scalar product of normal vector with velocity
179 for (unsigned i = 0; i < D; ++i) {
180 scalarVelocity += vectorVelocity[i] * normalVector[i];
181 }
182
183 if (scalarVelocity == T(0)) {
184 return {0, 0};
185 }
186
187 T hamiltonian =
188 Norm(calculateGradient(viennahrle::Index<D>(0))) * scalarVelocity;
189 T dissipation = 0.; // dissipation
190
191 // dissipation block
192 {
193 // reserve alphas for all points in local stencil
194 std::vector<VectorType<T, D>> alphas;
195 alphas.reserve(numStencilPoints);
196
197 viennahrle::Index<D> currentIndex(-order);
198 for (size_t i = 0; i < numStencilPoints; ++i) {
199 VectorType<T, D> alpha;
200 Vec3D<T> normal = calculateNormal(currentIndex);
201
202 // Check for corrupted normal
203 if ((std::abs(normal[0]) < 1e-6) && (std::abs(normal[1]) < 1e-6) &&
204 (std::abs(normal[2]) < 1e-6)) {
205 continue;
206 }
207
208 Vec3D<T> normal_p{normal[0], normal[1], normal[2]};
209 Vec3D<T> normal_n{normal[0], normal[1], normal[2]};
210
211 VectorType<T, D> velocityDelta;
212 std::fill(velocityDelta.begin(), velocityDelta.end(), 0.);
213
214 // get local velocity
215 Vec3D<T> localCoordArray = coordArray;
216 for (unsigned dir = 0; dir < D; ++dir)
217 localCoordArray[dir] += currentIndex[dir];
218
219 T localScalarVelocity = velocities->getScalarVelocity(
220 localCoordArray, material, normal_p,
221 neighborIterator.getCenter().getPointId());
222 Vec3D<T> localVectorVelocity = velocities->getVectorVelocity(
223 localCoordArray, material, normal_p,
224 neighborIterator.getCenter().getPointId());
225 // now calculate scalar product of normal vector with velocity
226 for (unsigned dir = 0; dir < D; ++dir) {
227 localScalarVelocity += localVectorVelocity[dir] * normal[dir];
228 }
229
230 const T DN = std::abs(normalEpsilon * scalarVelocity);
231
232 for (int k = 0; k < D; ++k) {
233
234 normal_p[k] -= DN; // p=previous
235 normal_n[k] += DN; // n==next
236
237 T vp = velocities->getScalarVelocity(
238 localCoordArray, material, normal_p,
239 neighborIterator.getCenter().getPointId());
240 T vn = velocities->getScalarVelocity(
241 localCoordArray, material, normal_n,
242 neighborIterator.getCenter().getPointId());
243 // central difference
244 velocityDelta[k] = (vn - vp) / (2.0 * DN);
245
246 normal_p[k] += DN;
247 normal_n[k] -= DN;
248 }
249
250 // determine \partial H / \partial phi_l
251 for (int k = 0; k < D; ++k) { // iterate over dimensions
252
253 // Monti term
254 T monti = 0;
255 // Toifl Quell term
256 T toifl = 0;
257
258 VectorType<T, D> gradient = calculateGradient(currentIndex);
259
260 for (int j = 0; j < D - 1; ++j) { // phi_p**2 + phi_q**2
261 int idx = (k + 1 + j) % D;
262 monti += gradient[idx] * gradient[idx];
263 toifl += gradient[idx] * velocityDelta[idx];
264 }
265 // denominator: |grad(phi)|^2
266 T denom = DotProduct(gradient, gradient);
267 monti *= velocityDelta[k] / denom;
268 toifl *= -gradient[k] / denom;
269
270 // Osher (constant V) term
271 T osher = localScalarVelocity * normal[k];
272 // Total derivative is sum of terms given above
273 alpha[k] = std::fabs(monti + toifl + osher);
274 }
275
276 alphas.push_back(alpha);
277
278 // increment current index
279 {
280 int dim = 0;
281 for (; dim < D - 1; ++dim) {
282 if (currentIndex[dim] < order)
283 break;
284 currentIndex[dim] = -order;
285 }
286 ++currentIndex[dim];
287 }
288 }
289
290 // determine max alphas for every axis
291 VectorType<T, D> gradientDiff = calculateGradientDiff();
292 for (int d = 0; d < D; ++d) {
293 T maxAlpha = 0;
294
295 for (size_t i = 0; i < numStencilPoints; ++i) {
296 maxAlpha = std::max(maxAlpha, alphas[i][d]);
297 }
298
299 finalAlphas[d] = maxAlpha;
300 dissipation += maxAlpha * gradientDiff[d];
301 }
302
303 return {hamiltonian, dissipation};
304 }
305 }
306
307 void reduceTimeStepHamiltonJacobi(double &MaxTimeStep,
308 double gridDelta) const {
309 constexpr double alpha_maxCFL = 1.0;
310 // second time step test, based on alphas
311
312 double timeStep = 0;
313 for (int i = 0; i < D; ++i) {
314 timeStep += finalAlphas[i] / gridDelta;
315 }
316
317 timeStep = alpha_maxCFL / timeStep;
318 MaxTimeStep = std::min(timeStep, MaxTimeStep);
319 }
320};
321} // namespace lsInternal
322
323namespace viennals {
324
325using namespace viennacore;
326
338template <class T, int D>
340 std::vector<SmartPointer<Domain<T, D>>> &levelSets,
341 std::vector<bool> isDepo) {
342 if (isDepo.size() < levelSets.size()) {
343 Logger::getInstance()
344 .addWarning(
345 "PrepareStencilLocalLaxFriedrichs: isDepo does not have enough "
346 "elements. Assuming all higher layers are not depo layers.")
347 .print();
348 }
349 // always resize, so it has the correct number of elements
350 isDepo.resize(levelSets.size(), false);
351
352 // Begin with the biggest level set (top LS wrapped around all others)
353 auto layerIt = levelSets.rbegin();
354
355 // make empty LS which will contain the final top layer
356 auto finalTop = SmartPointer<Domain<T, D>>::New(levelSets[0]->getGrid());
357
358 bool layerAboveIsDepo = false;
359 auto depoIt = isDepo.rbegin();
360
361 // subtract each first non-depo layer below a depo layer from the latter
362 // then union these results to the final layer
363 for (auto maskIt = levelSets.rbegin(); maskIt != levelSets.rend(); ++maskIt) {
364 if (isDepo[*depoIt]) {
365 // this layer will be deposited on
366 if (!layerAboveIsDepo) {
367 layerIt = maskIt;
368 layerAboveIsDepo = true;
369 }
370 } else {
371 // no depo on this layer
372 if (layerAboveIsDepo) {
373 auto layerAbove = SmartPointer<Domain<T, D>>::New(*layerIt);
374 BooleanOperation<T, D>(layerAbove, *maskIt,
376 .apply();
377
378 BooleanOperation<T, D>(finalTop, layerAbove,
380 .apply();
381 Prune<T, D>(finalTop).apply();
382 }
383 layerAboveIsDepo = false;
384 }
385 ++depoIt;
386 }
387
388 // If the lowest layer is a depo substrate, add it to the final top layer
389 if (layerAboveIsDepo) {
391 .apply();
392 }
393
395 levelSets.back()->deepCopy(finalTop);
396}
397
400template <class T, int D>
402 std::vector<SmartPointer<Domain<T, D>>> &levelSets) {
403 auto layerIt = levelSets.rbegin();
404 auto lastIt = ++levelSets.rbegin();
405
408 .apply();
409}
410
411} // namespace viennals
static constexpr unsigned getNumberOfValues()
Definition lsFiniteDifferences.hpp:40
static T calculateGradient(const T *values, const double &delta)
Definition lsFiniteDifferences.hpp:198
static T calculateGradientDiff(const T *values, const double &delta)
Definition lsFiniteDifferences.hpp:204
std::pair< T, T > operator()(const viennahrle::Index< D > &indices, int material)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:132
StencilLocalLaxFriedrichsScalar(LevelSetType passedlsDomain, SmartPointer< viennals::VelocityField< T > > vel, double a=1.0)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:122
static void prepareLS(LevelSetType passedlsDomain)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:116
void reduceTimeStepHamiltonJacobi(double &MaxTimeStep, double gridDelta) const
Definition lsStencilLocalLaxFriedrichsScalar.hpp:307
This class is used to perform boolean operations on two level sets and write the resulting level set ...
Definition lsBooleanOperation.hpp:43
void apply()
Perform operation.
Definition lsBooleanOperation.hpp:316
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
Expands the levelSet to the specified number of layers. The largest value in the levelset is thus wid...
Definition lsExpand.hpp:17
void apply()
Apply the expansion to the specified width.
Definition lsExpand.hpp:44
Removes all level set points, which do not have at least one oppositely signed neighbour (Meaning the...
Definition lsPrune.hpp:17
void apply()
removes all grid points, which do not have at least one opposite signed neighbour returns the number ...
Definition lsPrune.hpp:74
Abstract class defining the interface for the velocity field used during advection using lsAdvect.
Definition lsVelocityField.hpp:11
Definition lsCurvatureFormulas.hpp:9
DifferentiationSchemeEnum
Definition lsFiniteDifferences.hpp:10
@ FIRST_ORDER
Definition lsFiniteDifferences.hpp:11
Definition lsAdvect.hpp:36
void PrepareStencilLocalLaxFriedrichs(std::vector< SmartPointer< Domain< T, D > > > &levelSets, std::vector< bool > isDepo)
This function creates the specialized layer wrapping which produces better results for the SSLF integ...
Definition lsStencilLocalLaxFriedrichsScalar.hpp:339
void FinalizeStencilLocalLaxFriedrichs(std::vector< SmartPointer< Domain< T, D > > > &levelSets)
After advection using the SLLF layer wrapping approach is done, restore the original layer wrapping u...
Definition lsStencilLocalLaxFriedrichsScalar.hpp:401
@ INVERT
Definition lsBooleanOperation.hpp:29
@ RELATIVE_COMPLEMENT
Definition lsBooleanOperation.hpp:28
@ UNION
Definition lsBooleanOperation.hpp:27
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68