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 static double maxDissipation; // default: std::numeric_limits<double>::max();
44
45 static T pow2(const T &value) { return value * value; }
46
47 Vec3D<T> calculateNormal(const viennahrle::Index<D> &offset) {
48 Vec3D<T> normal = {0.0, 0.0, 0.0};
49 constexpr int startIndex = -1;
50 T modulus = 0.;
51
52 for (unsigned i = 0; i < D; ++i) {
53 viennahrle::Index<D> index(offset);
54 std::vector<T> values;
55 for (unsigned j = 0; j < 3; ++j) {
56 index[i] = startIndex + j;
57 values.push_back(neighborIterator.getNeighbor(index).getValue());
58 }
60 values.data(), levelSet->getGrid().getGridDelta());
61 modulus += normal[i] * normal[i];
62 }
63 modulus = std::sqrt(modulus);
64 for (unsigned i = 0; i < D; ++i) {
65 normal[i] /= modulus;
66 }
67 return normal;
68 }
69
70 VectorType<T, D> calculateGradient(const viennahrle::Index<D> &offset) {
71 VectorType<T, D> gradient;
72
73 constexpr unsigned numValues =
75 const int startIndex = -std::floor(numValues / 2);
76
77 for (unsigned i = 0; i < D; ++i) {
78 viennahrle::Index<D> index(offset);
79 std::vector<T> values;
80 for (unsigned j = 0; j < numValues; ++j) {
81 index[i] = startIndex + j;
82 values.push_back(neighborIterator.getNeighbor(index).getValue());
83 }
84
85 gradient[i] =
87 values.data(), levelSet->getGrid().getGridDelta());
88 }
89
90 return gradient;
91 }
92
93 VectorType<T, D> calculateGradientDiff() {
94 VectorType<T, D> gradient;
95
96 const unsigned numValues =
98 const int startIndex = -std::floor(numValues / 2);
99
100 for (unsigned i = 0; i < D; ++i) {
101 viennahrle::Index<D> index(0);
102 std::vector<T> values;
103 for (unsigned j = 0; j < numValues; ++j) {
104 index[i] = startIndex + j;
105 values.push_back(neighborIterator.getNeighbor(index).getValue());
106 }
107
108 gradient[i] =
110 &(values[0]), levelSet->getGrid().getGridDelta());
111 }
112
113 return gradient;
114 }
115
116public:
117 static void prepareLS(LevelSetType passedlsDomain) {
118 // Expansion of sparse field must depend on spatial derivative order
119 // AND slf stencil order! --> currently assume scheme = 3rd order always
120 viennals::Expand<T, D>(passedlsDomain, 2 * (order + 1) + 4).apply();
121 }
122
123 StencilLocalLaxFriedrichsScalar(LevelSetType passedlsDomain,
124 SmartPointer<viennals::VelocityField<T>> vel,
125 double a = 1.0)
126 : levelSet(passedlsDomain), velocities(vel),
127 neighborIterator(levelSet->getDomain()), alphaFactor(a) {
128 for (int i = 0; i < 3; ++i) {
129 finalAlphas[i] = 0;
130 }
131 }
132
133 static void setMaxDissipation(double maxDiss) { maxDissipation = maxDiss; }
134
135 std::pair<T, T> operator()(const viennahrle::Index<D> &indices,
136 int material) {
137 auto &grid = levelSet->getGrid();
138 double gridDelta = grid.getGridDelta();
139
140 Vec3D<T> coordinate{0., 0., 0.};
141 for (unsigned i = 0; i < D; ++i) {
142 coordinate[i] = indices[i] * gridDelta;
143 }
144
145 // move neighborIterator to current position
146 neighborIterator.goToIndicesSequential(indices);
147
148 // convert coordinate to std array for interface
149 Vec3D<T> coordArray{coordinate[0], coordinate[1], coordinate[2]};
150
151 // if there is a vector velocity, we need to project it onto a scalar
152 // velocity first using its normal vector
153 // /*if (vectorVelocity != Vec3D<T>({}))*/ {
154 Vec3D<T> normalVector;
155 T denominator = 0; // normal modulus
156 for (unsigned i = 0; i < D; i++) {
157 viennahrle::Index<D> neighborIndex(0);
158 neighborIndex[i] = 1;
159 // normal vector calculation
160 T pos = neighborIterator.getNeighbor(neighborIndex).getValue() -
161 neighborIterator.getCenter().getValue();
162 T neg = neighborIterator.getCenter().getValue() -
163 neighborIterator.getNeighbor(-neighborIndex).getValue();
164 normalVector[i] = (pos + neg) * 0.5;
165 // normalise normal vector
166 denominator += normalVector[i] * normalVector[i];
167 }
168 denominator = std::sqrt(denominator);
169
170 for (unsigned i = 0; i < D; ++i) {
171 normalVector[i] /= denominator;
172 }
173
174 double scalarVelocity = velocities->getScalarVelocity(
175 coordArray, material, normalVector,
176 neighborIterator.getCenter().getPointId());
177 auto vectorVelocity = velocities->getVectorVelocity(
178 coordArray, material, normalVector,
179 neighborIterator.getCenter().getPointId());
180
181 // now calculate scalar product of normal vector with velocity
182 for (unsigned i = 0; i < D; ++i) {
183 scalarVelocity += vectorVelocity[i] * normalVector[i];
184 }
185
186 if (scalarVelocity == T(0)) {
187 return {0, 0};
188 }
189
190 T hamiltonian =
191 Norm(calculateGradient(viennahrle::Index<D>(0))) * scalarVelocity;
192 T dissipation = 0.; // dissipation
193
194 // dissipation block
195 {
196 // reserve alphas for all points in local stencil
197 std::vector<VectorType<T, D>> alphas;
198 alphas.reserve(numStencilPoints);
199
200 viennahrle::Index<D> currentIndex(-order);
201 for (size_t i = 0; i < numStencilPoints; ++i) {
202 VectorType<T, D> alpha;
203 Vec3D<T> normal = calculateNormal(currentIndex);
204
205 // Check for corrupted normal
206 if ((std::abs(normal[0]) < 1e-6) && (std::abs(normal[1]) < 1e-6) &&
207 (std::abs(normal[2]) < 1e-6)) {
208 continue;
209 }
210
211 Vec3D<T> normal_p{normal[0], normal[1], normal[2]};
212 Vec3D<T> normal_n{normal[0], normal[1], normal[2]};
213
214 VectorType<T, D> velocityDelta;
215 std::fill(velocityDelta.begin(), velocityDelta.end(), 0.);
216
217 // get local velocity
218 Vec3D<T> localCoordArray = coordArray;
219 for (unsigned dir = 0; dir < D; ++dir)
220 localCoordArray[dir] += currentIndex[dir];
221
222 T localScalarVelocity = velocities->getScalarVelocity(
223 localCoordArray, material, normal_p,
224 neighborIterator.getCenter().getPointId());
225 Vec3D<T> localVectorVelocity = velocities->getVectorVelocity(
226 localCoordArray, material, normal_p,
227 neighborIterator.getCenter().getPointId());
228 // now calculate scalar product of normal vector with velocity
229 for (unsigned dir = 0; dir < D; ++dir) {
230 localScalarVelocity += localVectorVelocity[dir] * normal[dir];
231 }
232
233 const T DN = std::abs(normalEpsilon * scalarVelocity);
234
235 for (int k = 0; k < D; ++k) {
236
237 normal_p[k] -= DN; // p=previous
238 normal_n[k] += DN; // n==next
239
240 T vp = velocities->getScalarVelocity(
241 localCoordArray, material, normal_p,
242 neighborIterator.getCenter().getPointId());
243 T vn = velocities->getScalarVelocity(
244 localCoordArray, material, normal_n,
245 neighborIterator.getCenter().getPointId());
246 // central difference
247 velocityDelta[k] = (vn - vp) / (2.0 * DN);
248
249 normal_p[k] += DN;
250 normal_n[k] -= DN;
251 }
252
253 // determine \partial H / \partial phi_l
254 for (int k = 0; k < D; ++k) { // iterate over dimensions
255
256 // Monti term
257 T monti = 0;
258 // Toifl Quell term
259 T toifl = 0;
260
261 VectorType<T, D> gradient = calculateGradient(currentIndex);
262
263 for (int j = 0; j < D - 1; ++j) { // phi_p**2 + phi_q**2
264 int idx = (k + 1 + j) % D;
265 monti += gradient[idx] * gradient[idx];
266 toifl += gradient[idx] * velocityDelta[idx];
267 }
268 // denominator: |grad(phi)|^2
269 T denom = DotProduct(gradient, gradient);
270 monti *= velocityDelta[k] / denom;
271 toifl *= -gradient[k] / denom;
272
273 // Osher (constant V) term
274 T osher = localScalarVelocity * normal[k];
275 // Total derivative is sum of terms given above
276 alpha[k] = std::fabs(monti + toifl + osher);
277
278 if (alpha[k] > maxDissipation) {
279 alpha[k] = 0.;
280 }
281 }
282
283 alphas.push_back(alpha);
284
285 // increment current index
286 {
287 int dim = 0;
288 for (; dim < D - 1; ++dim) {
289 if (currentIndex[dim] < order)
290 break;
291 currentIndex[dim] = -order;
292 }
293 ++currentIndex[dim];
294 }
295 }
296
297 // determine max alphas for every axis
298 VectorType<T, D> gradientDiff = calculateGradientDiff();
299 for (int d = 0; d < D; ++d) {
300 T maxAlpha = 0;
301
302 for (size_t i = 0; i < numStencilPoints; ++i) {
303 maxAlpha = std::max(maxAlpha, alphas[i][d]);
304 }
305
306 finalAlphas[d] = maxAlpha;
307 dissipation += maxAlpha * gradientDiff[d];
308 }
309
310 return {hamiltonian, dissipation};
311 }
312 }
313
314 void reduceTimeStepHamiltonJacobi(double &MaxTimeStep,
315 double gridDelta) const {
316 constexpr double alpha_maxCFL = 1.0;
317 // second time step test, based on alphas
318
319 double timeStep = 0;
320 for (int i = 0; i < D; ++i) {
321 timeStep += finalAlphas[i] / gridDelta;
322 }
323
324 timeStep = alpha_maxCFL / timeStep;
325 MaxTimeStep = std::min(timeStep, MaxTimeStep);
326 }
327};
328
329template <class T, int D, int order,
330 DifferentiationSchemeEnum finiteDifferenceScheme>
331double StencilLocalLaxFriedrichsScalar<T, D, order,
332 finiteDifferenceScheme>::maxDissipation =
333 std::numeric_limits<double>::max();
334
335} // namespace lsInternal
336
337namespace viennals {
338
339using namespace viennacore;
340
352template <class T, int D>
354 std::vector<SmartPointer<Domain<T, D>>> &levelSets,
355 std::vector<bool> isDepo) {
356 if (isDepo.size() < levelSets.size()) {
357 Logger::getInstance()
358 .addWarning(
359 "PrepareStencilLocalLaxFriedrichs: isDepo does not have enough "
360 "elements. Assuming all higher layers are not depo layers.")
361 .print();
362 }
363 // always resize, so it has the correct number of elements
364 isDepo.resize(levelSets.size(), false);
365
366 // Begin with the biggest level set (top LS wrapped around all others)
367 auto layerIt = levelSets.rbegin();
368
369 // make empty LS which will contain the final top layer
370 auto finalTop = SmartPointer<Domain<T, D>>::New(levelSets[0]->getGrid());
371
372 bool layerAboveIsDepo = false;
373 auto depoIt = isDepo.rbegin();
374
375 // subtract each first non-depo layer below a depo layer from the latter
376 // then union these results to the final layer
377 for (auto maskIt = levelSets.rbegin(); maskIt != levelSets.rend(); ++maskIt) {
378 if (isDepo[*depoIt]) {
379 // this layer will be deposited on
380 if (!layerAboveIsDepo) {
381 layerIt = maskIt;
382 layerAboveIsDepo = true;
383 }
384 } else {
385 // no depo on this layer
386 if (layerAboveIsDepo) {
387 auto layerAbove = SmartPointer<Domain<T, D>>::New(*layerIt);
388 BooleanOperation<T, D>(layerAbove, *maskIt,
390 .apply();
391
392 BooleanOperation<T, D>(finalTop, layerAbove,
394 .apply();
395 Prune<T, D>(finalTop).apply();
396 }
397 layerAboveIsDepo = false;
398 }
399 ++depoIt;
400 }
401
402 // If the lowest layer is a depo substrate, add it to the final top layer
403 if (layerAboveIsDepo) {
405 .apply();
406 }
407
409 levelSets.back()->deepCopy(finalTop);
410}
411
414template <class T, int D>
416 std::vector<SmartPointer<Domain<T, D>>> &levelSets) {
417 auto layerIt = levelSets.rbegin();
418 auto lastIt = ++levelSets.rbegin();
419
422 .apply();
423}
424
425} // namespace viennals
constexpr int D
Definition Epitaxy.cpp:9
double T
Definition Epitaxy.cpp:10
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:133
std::pair< T, T > operator()(const viennahrle::Index< D > &indices, int material)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:135
StencilLocalLaxFriedrichsScalar(LevelSetType passedlsDomain, SmartPointer< viennals::VelocityField< T > > vel, double a=1.0)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:123
static void prepareLS(LevelSetType passedlsDomain)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:117
void reduceTimeStepHamiltonJacobi(double &MaxTimeStep, double gridDelta) const
Definition lsStencilLocalLaxFriedrichsScalar.hpp:314
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:332
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:353
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:415
@ INVERT
Definition lsBooleanOperation.hpp:29
@ RELATIVE_COMPLEMENT
Definition lsBooleanOperation.hpp:28
@ UNION
Definition lsBooleanOperation.hpp:27