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
12// Choose adaptive epsilon strategy
13// 0: Original fixed epsilon
14// 1: Simple adaptive epsilon (recommended)
15// 2: Full adaptive epsilon (most accurate but slower)
16#ifndef ADAPTIVE_EPSILON_STRATEGY
17#define ADAPTIVE_EPSILON_STRATEGY 0
18#endif
19
20namespace lsInternal {
21
22using namespace viennacore;
23
30template <class T, int D, int order,
31 DifferentiationSchemeEnum finiteDifferenceScheme =
34 using LevelSetType = SmartPointer<viennals::Domain<T, D>>;
35 using LevelSetsType = std::vector<LevelSetType>;
36
37 LevelSetType levelSet;
38 SmartPointer<viennals::VelocityField<T>> velocities;
39 viennahrle::SparseBoxIterator<viennahrle::Domain<T, D>,
40 static_cast<int>(finiteDifferenceScheme) + 1 +
41 order>
42 neighborIterator;
43 const double alphaFactor;
44 const double baseNormalEpsilon =
45 std::cbrt(std::numeric_limits<double>::epsilon());
46
47 // Final dissipation coefficients that are used by the time integrator. If
48 // D==2 last entries are 0.
49 Vec3D<T> finalAlphas;
50 static constexpr unsigned numStencilPoints = hrleUtil::pow(2 * order + 1, D);
51 static double maxDissipation; // default: std::numeric_limits<double>::max();
52
53 // Adaptive epsilon control parameters
54 static constexpr T velocityScaleFactor =
55 0.5; // Controls velocity-based scaling
56 static constexpr T smoothnessThreshold =
57 0.1; // Threshold for smoothness detection
58 static constexpr T epsilonMinScale = 0.1; // Minimum epsilon scale factor
59 static constexpr T epsilonMaxScale = 100.0; // Maximum epsilon scale factor
60
61 static T pow2(const T &value) { return value * value; }
62
63 // Calculate adaptive epsilon based on velocity magnitude, grid resolution,
64 // and local velocity smoothness
65 T calculateAdaptiveEpsilon(T velocityMagnitude, T gridDelta,
66 const Vec3D<T> &localCoordArray, int material,
67 const Vec3D<T> &normal) const {
68 // Base epsilon scaled by grid resolution
69 T adaptiveEpsilon = baseNormalEpsilon * gridDelta;
70
71 // Scale with velocity magnitude to maintain numerical stability
72 // Use larger epsilon for larger velocities to avoid catastrophic
73 // cancellation
74 T absVelocity = std::abs(velocityMagnitude);
75 if (absVelocity > 1e-12) {
76 T velocityScale = 1.0 + velocityScaleFactor * std::log(1.0 + absVelocity);
77 adaptiveEpsilon *= velocityScale;
78 }
79
80 // Estimate local velocity smoothness using a simplified approach
81 // Only check smoothness if velocity is significant
82 if (absVelocity > 1e-10) {
83 T maxRelativeVariation = 0.0;
84
85 // Sample velocity in each coordinate direction with small perturbation
86 for (int dir = 0; dir < D; ++dir) {
87 Vec3D<T> perturbedNormal = normal;
88 perturbedNormal[dir] += smoothnessThreshold;
89
90 // Normalize the perturbed normal
91 T normSq = 0.0;
92 for (int i = 0; i < D; ++i) {
93 normSq += perturbedNormal[i] * perturbedNormal[i];
94 }
95 if (normSq > 1e-12) {
96 T invNorm = 1.0 / std::sqrt(normSq);
97 for (int i = 0; i < D; ++i) {
98 perturbedNormal[i] *= invNorm;
99 }
100
101 T perturbedVel = velocities->getScalarVelocity(
102 localCoordArray, material, perturbedNormal,
103 neighborIterator.getCenter().getPointId());
104
105 T relativeVariation =
106 std::abs(perturbedVel - velocityMagnitude) / absVelocity;
107 maxRelativeVariation =
108 std::max(maxRelativeVariation, relativeVariation);
109 }
110 }
111
112 // If velocity varies significantly, use larger epsilon for stability
113 if (maxRelativeVariation > smoothnessThreshold) {
114 T smoothnessFactor = 1.0 + maxRelativeVariation;
115 adaptiveEpsilon *= smoothnessFactor;
116 }
117 }
118
119 // Clamp epsilon to reasonable bounds
120 T minEpsilon = baseNormalEpsilon * gridDelta * epsilonMinScale;
121 T maxEpsilon = baseNormalEpsilon * gridDelta * epsilonMaxScale;
122
123 return std::max(minEpsilon, std::min(maxEpsilon, adaptiveEpsilon));
124 }
125
126 Vec3D<T> calculateNormal(const viennahrle::Index<D> &offset) {
127 Vec3D<T> normal = {0.0, 0.0, 0.0};
128 constexpr int startIndex = -1;
129 T modulus = 0.;
130
131 for (unsigned i = 0; i < D; ++i) {
132 viennahrle::Index<D> index(offset);
133 std::array<T, 3> values;
134 for (unsigned j = 0; j < 3; ++j) {
135 index[i] = startIndex + j;
136 values[j] = neighborIterator.getNeighbor(index).getValue();
137 }
139 values.data(), levelSet->getGrid().getGridDelta());
140 modulus += normal[i] * normal[i];
141 }
142 modulus = 1 / std::sqrt(modulus);
143 for (unsigned i = 0; i < D; ++i) {
144 normal[i] *= modulus;
145 }
146 return normal;
147 }
148
149 VectorType<T, D> calculateGradient(const viennahrle::Index<D> &offset) {
150 VectorType<T, D> gradient;
151
152 constexpr unsigned numValues =
154 const int startIndex = -std::floor(numValues / 2);
155
156 for (unsigned i = 0; i < D; ++i) {
157 viennahrle::Index<D> index(offset);
158 std::array<T, numValues> values;
159 for (unsigned j = 0; j < numValues; ++j) {
160 index[i] = startIndex + j;
161 values[j] = neighborIterator.getNeighbor(index).getValue();
162 }
163
164 gradient[i] =
166 values.data(), levelSet->getGrid().getGridDelta());
167 }
168
169 return gradient;
170 }
171
172 VectorType<T, D> calculateGradientDiff() {
173 VectorType<T, D> gradient;
174
175 constexpr unsigned numValues =
177 const int startIndex =
178 -std::floor(numValues / 2); // std::floor constexpr in C++23
179
180 for (unsigned i = 0; i < D; ++i) {
181 viennahrle::Index<D> index(0);
182 std::array<T, numValues> values;
183 for (unsigned j = 0; j < numValues; ++j) {
184 index[i] = startIndex + j;
185 values[j] = neighborIterator.getNeighbor(index).getValue();
186 }
187
188 gradient[i] =
190 values.data(), levelSet->getGrid().getGridDelta());
191 }
192
193 return gradient;
194 }
195
196public:
197 static void prepareLS(LevelSetType passedlsDomain) {
198 // Expansion of sparse field must depend on spatial derivative order
199 // AND slf stencil order! --> currently assume scheme = 3rd order always
200 viennals::Expand<T, D>(passedlsDomain, 2 * (order + 1) + 4).apply();
201 }
202
203 StencilLocalLaxFriedrichsScalar(LevelSetType passedlsDomain,
204 SmartPointer<viennals::VelocityField<T>> vel,
205 double a = 1.0)
206 : levelSet(passedlsDomain), velocities(vel),
207 neighborIterator(levelSet->getDomain()), alphaFactor(a) {
208 for (int i = 0; i < 3; ++i) {
209 finalAlphas[i] = 0;
210 }
211 }
212
213 static void setMaxDissipation(double maxDiss) { maxDissipation = maxDiss; }
214
215 // Simple adaptive epsilon calculation for better performance
216 T calculateSimpleAdaptiveEpsilon(T velocityMagnitude, T gridDelta) const {
217 T absVel = std::abs(velocityMagnitude);
218 T baseEps = baseNormalEpsilon * gridDelta;
219
220 // Simple velocity-based scaling to avoid numerical issues
221 if (absVel > 1e-10) {
222 // Scale epsilon with velocity magnitude (clamped to reasonable range)
223 T scale = std::max(T(0.1), std::min(T(10.0), absVel));
224 return baseEps * scale;
225 }
226
227 return baseEps;
228 }
229
230 std::pair<T, T> operator()(const viennahrle::Index<D> &indices,
231 int material) {
232 auto &grid = levelSet->getGrid();
233 double gridDelta = grid.getGridDelta();
234
235 Vec3D<T> coordinate{0., 0., 0.};
236 for (unsigned i = 0; i < D; ++i) {
237 coordinate[i] = indices[i] * gridDelta;
238 }
239
240 // move neighborIterator to current position
241 neighborIterator.goToIndicesSequential(indices);
242
243 // if there is a vector velocity, we need to project it onto a scalar
244 // velocity first using its normal vector
245 // /*if (vectorVelocity != Vec3D<T>({}))*/ {
246 Vec3D<T> normalVector;
247 T denominator = 0; // normal modulus
248 for (unsigned i = 0; i < D; i++) {
249 viennahrle::Index<D> neighborIndex(0);
250 neighborIndex[i] = 1;
251 // normal vector calculation
252 T pos = neighborIterator.getNeighbor(neighborIndex).getValue() -
253 neighborIterator.getCenter().getValue();
254 T neg = neighborIterator.getCenter().getValue() -
255 neighborIterator.getNeighbor(-neighborIndex).getValue();
256 normalVector[i] = (pos + neg) * 0.5;
257 // normalise normal vector
258 denominator += normalVector[i] * normalVector[i];
259 }
260 denominator = 1 / std::sqrt(denominator);
261 for (unsigned i = 0; i < D; ++i) {
262 normalVector[i] *= denominator;
263 }
264
265 double scalarVelocity = velocities->getScalarVelocity(
266 coordinate, material, normalVector,
267 neighborIterator.getCenter().getPointId());
268 auto vectorVelocity = velocities->getVectorVelocity(
269 coordinate, material, normalVector,
270 neighborIterator.getCenter().getPointId());
271
272 // now calculate scalar product of normal vector with velocity
273 for (unsigned i = 0; i < D; ++i) {
274 scalarVelocity += vectorVelocity[i] * normalVector[i];
275 }
276
277 if (scalarVelocity == T(0)) {
278 return {0, 0};
279 }
280
281 T hamiltonian =
282 Norm(calculateGradient(viennahrle::Index<D>(0))) * scalarVelocity;
283 T dissipation = 0.; // dissipation
284
285 // dissipation block
286 {
287 // reserve alphas for all points in local stencil
288 std::array<VectorType<T, D>, numStencilPoints> alphas{};
289
290 viennahrle::Index<D> currentIndex(-order);
291 for (size_t i = 0; i < numStencilPoints; ++i) {
292 VectorType<T, D> alpha;
293 Vec3D<T> normal = calculateNormal(currentIndex);
294
295 // Check for corrupted normal
296 if ((std::abs(normal[0]) < 1e-6) && (std::abs(normal[1]) < 1e-6) &&
297 (std::abs(normal[2]) < 1e-6)) {
298 alphas[i].fill(0);
299 continue;
300 }
301
302 Vec3D<T> normal_p{normal[0], normal[1], normal[2]};
303 Vec3D<T> normal_n{normal[0], normal[1], normal[2]};
304
305 VectorType<T, D> velocityDelta;
306 velocityDelta.fill(0.);
307
308 // get local velocity
309 Vec3D<T> localCoordArray = coordinate;
310 for (unsigned dir = 0; dir < D; ++dir)
311 localCoordArray[dir] += currentIndex[dir];
312
313 T localScalarVelocity = velocities->getScalarVelocity(
314 localCoordArray, material, normal_p,
315 neighborIterator.getCenter().getPointId());
316 Vec3D<T> localVectorVelocity = velocities->getVectorVelocity(
317 localCoordArray, material, normal_p,
318 neighborIterator.getCenter().getPointId());
319 // now calculate scalar product of normal vector with velocity
320 for (unsigned dir = 0; dir < D; ++dir) {
321 localScalarVelocity += localVectorVelocity[dir] * normal[dir];
322 }
323
324 // Calculate epsilon based on selected strategy
325 T DN;
326#if ADAPTIVE_EPSILON_STRATEGY == 0
327 // Original fixed epsilon approach
328 DN = std::abs(baseNormalEpsilon * scalarVelocity);
329#elif ADAPTIVE_EPSILON_STRATEGY == 1
330 // Simple adaptive epsilon (recommended)
331 DN = calculateSimpleAdaptiveEpsilon(localScalarVelocity,
332 levelSet->getGrid().getGridDelta());
333#elif ADAPTIVE_EPSILON_STRATEGY == 2
334 // Full adaptive epsilon (most accurate but slower)
335 DN = calculateAdaptiveEpsilon(localScalarVelocity,
336 levelSet->getGrid().getGridDelta(),
337 localCoordArray, material, normal);
338#else
339 // Fallback to original method
340 DN = std::abs(baseNormalEpsilon * scalarVelocity);
341#endif
342
343 for (int k = 0; k < D; ++k) {
344
345 normal_p[k] -= DN; // p=previous
346 normal_n[k] += DN; // n=next
347
348 T vp = velocities->getScalarVelocity(
349 localCoordArray, material, normal_p,
350 neighborIterator.getCenter().getPointId());
351 T vn = velocities->getScalarVelocity(
352 localCoordArray, material, normal_n,
353 neighborIterator.getCenter().getPointId());
354 // central difference
355 velocityDelta[k] = (vn - vp) / (2.0 * DN);
356
357 normal_p[k] += DN;
358 normal_n[k] -= DN;
359 }
360
361 // determine \partial H / \partial phi_l
362 for (int k = 0; k < D; ++k) { // iterate over dimensions
363
364 // Monti term
365 T monti = 0;
366 // Toifl Quell term
367 T toifl = 0;
368
369 VectorType<T, D> gradient = calculateGradient(currentIndex);
370
371 for (int j = 0; j < D - 1; ++j) { // phi_p**2 + phi_q**2
372 int idx = (k + 1 + j) % D;
373 monti += gradient[idx] * gradient[idx];
374 toifl += gradient[idx] * velocityDelta[idx];
375 }
376 // denominator: |grad(phi)|^2
377 T denom = DotProduct(gradient, gradient);
378 monti *= velocityDelta[k] / denom;
379 toifl *= -gradient[k] / denom;
380
381 // Osher (constant V) term
382 T osher = localScalarVelocity * normal[k];
383 // Total derivative is sum of terms given above
384 alpha[k] = std::fabs(monti + toifl + osher);
385
386 if (alpha[k] > maxDissipation) {
387 alpha[k] = maxDissipation;
388 }
389 }
390
391 alphas[i] = alpha;
392
393 // increment current index
394 {
395 int dim = 0;
396 for (; dim < D - 1; ++dim) {
397 if (currentIndex[dim] < order)
398 break;
399 currentIndex[dim] = -order;
400 }
401 ++currentIndex[dim];
402 }
403 }
404
405 // determine max alphas for every axis
406 VectorType<T, D> gradientDiff = calculateGradientDiff();
407 for (int d = 0; d < D; ++d) {
408 T maxAlpha = 0;
409
410 for (size_t i = 0; i < numStencilPoints; ++i) {
411 maxAlpha = std::max(maxAlpha, alphas[i][d]);
412 }
413
414 finalAlphas[d] = maxAlpha;
415 dissipation += maxAlpha * gradientDiff[d];
416 }
417
418 return {hamiltonian, dissipation};
419 }
420 }
421
422 void reduceTimeStepHamiltonJacobi(double &MaxTimeStep,
423 double gridDelta) const {
424 constexpr double alpha_maxCFL = 1.0;
425 // second time step test, based on alphas
426
427 double timeStep = 0;
428 for (int i = 0; i < D; ++i) {
429 timeStep += finalAlphas[i] / gridDelta;
430 }
431
432 timeStep = alpha_maxCFL / timeStep;
433 MaxTimeStep = std::min(timeStep, MaxTimeStep);
434 }
435};
436
437template <class T, int D, int order,
438 DifferentiationSchemeEnum finiteDifferenceScheme>
439double StencilLocalLaxFriedrichsScalar<T, D, order,
440 finiteDifferenceScheme>::maxDissipation =
441 std::numeric_limits<double>::max();
442
443} // namespace lsInternal
444
445namespace viennals {
446
447using namespace viennacore;
448
460template <class T, int D>
462 std::vector<SmartPointer<Domain<T, D>>> &levelSets,
463 std::vector<bool> isDepo) {
464 if (isDepo.size() < levelSets.size()) {
465 Logger::getInstance()
466 .addWarning(
467 "PrepareStencilLocalLaxFriedrichs: isDepo does not have enough "
468 "elements. Assuming all higher layers are not depo layers.")
469 .print();
470 }
471 // always resize, so it has the correct number of elements
472 isDepo.resize(levelSets.size(), false);
473
474 // Begin with the biggest level set (top LS wrapped around all others)
475 auto layerIt = levelSets.rbegin();
476
477 // make empty LS which will contain the final top layer
478 auto finalTop = SmartPointer<Domain<T, D>>::New(levelSets[0]->getGrid());
479
480 bool layerAboveIsDepo = false;
481 auto depoIt = isDepo.rbegin();
482
483 // subtract each first non-depo layer below a depo layer from the latter
484 // then union these results to the final layer
485 for (auto maskIt = levelSets.rbegin(); maskIt != levelSets.rend(); ++maskIt) {
486 if (isDepo[*depoIt]) {
487 // this layer will be deposited on
488 if (!layerAboveIsDepo) {
489 layerIt = maskIt;
490 layerAboveIsDepo = true;
491 }
492 } else {
493 // no depo on this layer
494 if (layerAboveIsDepo) {
495 auto layerAbove = SmartPointer<Domain<T, D>>::New(*layerIt);
496 BooleanOperation<T, D>(layerAbove, *maskIt,
498 .apply();
499
500 BooleanOperation<T, D>(finalTop, layerAbove,
502 .apply();
503 Prune<T, D>(finalTop).apply();
504 }
505 layerAboveIsDepo = false;
506 }
507 ++depoIt;
508 }
509
510 // If the lowest layer is a depo substrate, add it to the final top layer
511 if (layerAboveIsDepo) {
513 .apply();
514 }
515
517 levelSets.back()->deepCopy(finalTop);
518}
519
522template <class T, int D>
524 std::vector<SmartPointer<Domain<T, D>>> &levelSets) {
525 auto layerIt = levelSets.rbegin();
526 auto lastIt = ++levelSets.rbegin();
527
530 .apply();
531}
532
533} // namespace viennals
constexpr int D
Definition Epitaxy.cpp:11
double T
Definition Epitaxy.cpp:12
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
static void setMaxDissipation(double maxDiss)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:213
std::pair< T, T > operator()(const viennahrle::Index< D > &indices, int material)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:230
StencilLocalLaxFriedrichsScalar(LevelSetType passedlsDomain, SmartPointer< viennals::VelocityField< T > > vel, double a=1.0)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:203
T calculateSimpleAdaptiveEpsilon(T velocityMagnitude, T gridDelta) const
Definition lsStencilLocalLaxFriedrichsScalar.hpp:216
static void prepareLS(LevelSetType passedlsDomain)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:197
void reduceTimeStepHamiltonJacobi(double &MaxTimeStep, double gridDelta) const
Definition lsStencilLocalLaxFriedrichsScalar.hpp:422
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
double StencilLocalLaxFriedrichsScalar< T, D, order, finiteDifferenceScheme >::maxDissipation
Definition lsStencilLocalLaxFriedrichsScalar.hpp:440
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:461
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:523
@ INVERT
Definition lsBooleanOperation.hpp:29
@ RELATIVE_COMPLEMENT
Definition lsBooleanOperation.hpp:28
@ UNION
Definition lsBooleanOperation.hpp:27