Loading [MathJax]/extensions/tex2jax.js
MAIA bb96820c
Multiphysics at AIA
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
lbinterfacedxqy.cpp
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#include "lbinterfacedxqy.h"
8
9#include <algorithm>
10#include "IO/infoout.h"
11#include "MEMORY/collector.h"
12#include "UTIL/parallelfor.h"
13#include "lbconstants.h"
14#include "lbfunctions.h"
15#include "lbinterfacecell.h"
16#include "lbparentcell.h"
17#include "lbsolver.h"
18
19using namespace std;
20using namespace lbconstants;
21
29template <MInt nDim, MInt nDist, class SysEqn>
30LbInterfaceDxQy<nDim, nDist, SysEqn>::LbInterfaceDxQy(LbSolver<nDim>* solver)
31 : LbInterface<nDim>(solver)
32
33{
34 m_solver = static_cast<LbSolverDxQy<nDim, nDist, SysEqn>*>(solver);
35
52 const MInt solverId = m_solver->solverId();
53 m_interfaceMethod = "FILIPPOVA";
54 m_interfaceMethod = Context::getSolverProperty<MString>("interfaceMethod", solverId, AT_, &m_interfaceMethod);
55
56 setInterfaceFunctions();
57
58 m_adaptationInitMethod = "INIT_DUPUIS_FILIPPOVA";
59 m_adaptationInitMethod =
60
61 Context::getSolverProperty<MString>("interfaceMethod", solverId, AT_, &m_adaptationInitMethod);
62 setAdaptationFunctions();
63}
64
69template <MInt nDim, MInt nDist, class SysEqn>
70LbInterfaceDxQy<nDim, nDist, SysEqn>::~LbInterfaceDxQy() = default;
71
77template <MInt nDim, MInt nDist, class SysEqn>
78void LbInterfaceDxQy<nDim, nDist, SysEqn>::setInterfaceFunctions() {
79 TRACE();
80
81 switch(string2enum(m_interfaceMethod)) {
82 // Default property value
83 case FILIPPOVA: {
84 if(m_isThermal) {
85 fProlongation = &LbInterfaceDxQy::prolongationThermalDupuis;
86 fRestriction = &LbInterfaceDxQy::restrictionThermalDupuis;
87 } else if(m_solver->isCompressible()) {
88 fProlongation = &LbInterfaceDxQy::prolongationDupuisCompressible;
89 fRestriction = &LbInterfaceDxQy::restrictionDupuisCompressible;
90 } else {
91 fProlongation = &LbInterfaceDxQy::prolongationDupuis;
92 fRestriction = &LbInterfaceDxQy::restrictionDupuis;
93 }
94 break;
95 }
96 case ROHDE: {
97 if(m_isThermal) {
98 fProlongation = &LbInterfaceDxQy::prolongationThermalRohde;
99 fRestriction = &LbInterfaceDxQy::restrictionThermalRohde;
100 } else {
101 fProlongation = &LbInterfaceDxQy::prolongationRohde;
102 fRestriction = &LbInterfaceDxQy::restrictionRohde;
103 }
104 break;
105 }
106
107 default: {
108 stringstream errorMessage;
109 errorMessage << " LbInterface::setInterfaceFunctions: Specified interface condition: " << m_interfaceMethod
110 << " does not exist. Exiting!";
111 mTerm(1, AT_, errorMessage.str());
112 }
113 }
114}
115
121template <MInt nDim, MInt nDist, class SysEqn>
122void LbInterfaceDxQy<nDim, nDist, SysEqn>::prolongation() {
123 (this->*fProlongation)();
124}
125
131template <MInt nDim, MInt nDist, class SysEqn>
132void LbInterfaceDxQy<nDim, nDist, SysEqn>::restriction() {
133 (this->*fRestriction)();
134}
135
142template <MInt nDim, MInt nDist, class SysEqn>
143void LbInterfaceDxQy<nDim, nDist, SysEqn>::prolongation10() {
144 MFloat& tmp2 = m_static_prolongation10_tmp2;
145 MFloat(&b)[2 * nDim] = m_static_prolongation10_b;
146 MFloat(&c)[nDim * nDim] = m_static_prolongation10_c;
147 MFloat& trace = m_static_prolongation10_trace;
148
149 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
150 level++) { // minLevel+level is the actual coarse resolution
151 // minLevel+level+1 is the actual fine resolution
152
153 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level)
154 == 0) { // perform prolongation WITHOUT temporal interpolation in every second timestep on finest level,
155 // in every fourth timestep on the next coarser level, and so on; starting from t = 1
156
157 for(MInt id = 0; id < m_interfaceChildren[level]->size(); id++) {
158 const MInt currentId = m_interfaceChildren[level]->a[id].m_cellId;
159 if(currentId >= m_solver->noInternalCells()) continue;
160 const MInt parentId = m_solver->c_parentId(m_interfaceChildren[level]->a[id].m_cellId);
161 if(parentId >= m_solver->noInternalCells()) continue;
162
163 const MFloat omega_F =
164 2.0
165 / (1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
166 // omega_C = 2.0 / ( 1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel()
167 // - level));
168
169 //--------------------------------------------------------------------------------------------
170 // Spatial interpolation and temporal extrapolation of macroscopic variables and stress tensor
171 // f(t+1) = 2*f(t) - f(t-1)
172
173 MFloat rho = 0.0;
174 std::array<MFloat, nDim> u{};
175 for(MInt m = 0; m < IPOW2(nDim); m++) {
176 rho += 2.0 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
177 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
178 u[0] += 2.0 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
179 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->U);
180 u[1] += 2.0 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
181 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->V);
182 IF_CONSTEXPR(nDim == 3)
183 u[2] += 2.0 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
184 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->W);
185 }
186
187 for(MInt m = 0; m < IPOW2(nDim); m++) {
188 rho -= m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
189 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
190 u[0] -= m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
191 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->U);
192 u[1] -= m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
193 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->V);
194 IF_CONSTEXPR(nDim == 3)
195 u[2] -= m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
196 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->W);
197 }
198
199 // compute strain rate tensor on fine grid: (du/dx)_fine = 1/2*(du/dx)_coarse
200 IF_CONSTEXPR(nDim == 3) {
201 for(MInt j = 0; j < 3; j++) {
202 // du[j]/dx
203 c[j + 0 * 3] =
204 2.0 * rho * F1B2
205 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0]
206 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
207 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j)
208 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
209 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2]
210 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3])
211 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)
212 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j))
213 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4]
214 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5])
215 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)
216 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j))
217 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
218 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7])
219 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
220 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)));
221 // du[j]/dy
222 c[j + 1 * 3] =
223 2.0 * rho * F1B2
224 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2]
225 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0])
226 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j)
227 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
228 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3]
229 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
230 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)
231 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j))
232 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
233 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4])
234 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)
235 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j))
236 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7]
237 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5])
238 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
239 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)));
240 // du[j]/dz
241 c[j + 2 * 3] =
242 2.0 * rho * F1B2
243 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4]
244 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0])
245 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j)
246 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
247 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5]
248 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
249 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)
250 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j))
251 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
252 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2])
253 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)
254 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j))
255 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7]
256 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3])
257 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
258 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)));
259 }
260
261 for(MInt j = 0; j < 3; j++) {
262 // du[j]/dx
263 c[j + 0 * 3] -=
264 1.0 * rho * F1B2
265 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0]
266 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
267 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j)
268 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
269 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2]
270 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3])
271 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)
272 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j))
273 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4]
274 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5])
275 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)
276 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j))
277 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
278 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7])
279 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
280 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6],
281 j)));
282 // du[j]/dy
283 c[j + 1 * 3] -=
284 1.0 * rho * F1B2
285 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2]
286 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0])
287 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j)
288 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
289 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3]
290 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
291 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)
292 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j))
293 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
294 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4])
295 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)
296 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j))
297 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7]
298 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5])
299 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
300 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5],
301 j)));
302 // du[j]/dz
303 c[j + 2 * 3] -=
304 1.0 * rho * F1B2
305 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4]
306 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0])
307 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j)
308 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
309 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5]
310 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
311 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)
312 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j))
313 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
314 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2])
315 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)
316 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j))
317 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7]
318 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3])
319 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
320 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3],
321 j)));
322 }
323 trace = c[0] + c[4] + c[8];
324 }
325 else {
326 cout << "Strain rate tensor not implemented for 2D yet" << endl;
327 trace = F0;
328 }
329
330 //--------------------------------------------------------------------------------------------
331
332
333 // m_solver->a_variable(currentId, PV->RHO) = rho;
334 // m_solver->a_variable(currentId, PV->U) = u[0];
335 // m_solver->a_variable(currentId, PV->V) = u[1];
336 // m_solver->a_variable(currentId, PV->W) = u[2];
337
338 tmp2 = F0;
339 for(MInt d = 0; d < nDim; d++) {
340 tmp2 += (u[d] * u[d]);
341 b[2 * d] = -u[d];
342 b[2 * d + 1] = u[d];
343 }
344
345 // Calculate all equilibrium distributions
346 std::array<MFloat, nDist> eqDist;
347 eqDist = m_solver->getEqDists(rho, tmp2, u.data());
348
349 // interpolate and transform
350 for(MInt dist = 0; dist < nDist - 1; dist++) {
351 // only the missing distributions
352 if(m_solver->a_hasNeighbor(currentId, Ld::oppositeDist(dist)) == 0) {
353 const MFloat tp = Ld::tp(Ld::distType(dist));
354
355 // put together eq-part and rescaled non-eq part; scale factor = (omega_C/omega_F) * F1B2
356
357 m_solver->a_oldDistribution(currentId, dist) = eqDist[dist] + tp * trace / omega_F;
358
359 for(MInt k = 0; k < nDim; k++) {
360 for(MInt l = 0; l < nDim; l++) {
361 m_solver->a_oldDistribution(currentId, dist) -=
362 tp * F1BCSsq * (Ld::idFld(dist, k) - 1) * (Ld::idFld(dist, l) - 1) * c[l + nDim * k] / omega_F;
363 }
364 }
365 }
366 }
367 } // end of loop over all interface cells
368
369 }
370 // If NOT in sync with parent cells apply linear time interpolation
371 else {
372 if((globalTimeStep - 1) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level)
373 == 0) { // perform prolongation WITH temporal interpolation in every second timestep on finest level,
374 // in every fourth timestep on the next coarser level, and so on; starting from t = 2
375
376 for(MInt id = 0; id < m_interfaceChildren[level]->size(); id++) {
377 const MInt currentId = m_interfaceChildren[level]->a[id].m_cellId;
378 if(currentId >= m_solver->noInternalCells()) continue;
379 const MInt parentId = m_solver->c_parentId(m_interfaceChildren[level]->a[id].m_cellId);
380 if(parentId >= m_solver->noInternalCells()) continue;
381
382 const MFloat omega_F =
383 2.0
384 / (1.0
385 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
386
387 //--------------------------------------------------------------------------------------------
388 // Spatial interpolation and temporal extrapolation of macroscopic variables and stress tensor
389 // f(t+1/2) = 1.5*f(t) - 0.5*f(t-1)
390
391 MFloat rho = 0;
392 std::array<MFloat, nDim> u{};
393
394 for(MInt m = 0; m < IPOW2(nDim); m++) {
395 rho += 1.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
396 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
397 u[0] += 1.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
398 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->U);
399 u[1] += 1.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
400 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->V);
401 IF_CONSTEXPR(nDim == 3)
402 u[2] += 1.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
403 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->W);
404 }
405
406 for(MInt m = 0; m < IPOW2(nDim); m++) {
407 rho -= 0.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
408 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
409 u[0] -= 0.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
410 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->U);
411 u[1] -= 0.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
412 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->V);
413 IF_CONSTEXPR(nDim == 3)
414 u[2] -= 0.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
415 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->W);
416 }
417
418 // compute strain rate tensor on fine grid: (du/dx)_fine = 1/2*(du/dx)_coarse
419 IF_CONSTEXPR(nDim == 3) {
420 for(MInt j = 0; j < 3; j++) {
421 // du[j]/dx
422 c[j + 0 * 3] =
423 1.5 * rho * F1B2
424 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0]
425 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
426 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j)
427 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
428 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2]
429 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3])
430 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)
431 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j))
432 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4]
433 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5])
434 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)
435 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j))
436 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
437 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7])
438 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
439 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6],
440 j)));
441 // du[j]/dy
442 c[j + 1 * 3] =
443 1.5 * rho * F1B2
444 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2]
445 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0])
446 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j)
447 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
448 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3]
449 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
450 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)
451 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j))
452 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
453 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4])
454 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)
455 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j))
456 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7]
457 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5])
458 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
459 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5],
460 j)));
461 // du[j]/dz
462 c[j + 2 * 3] =
463 1.5 * rho * F1B2
464 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4]
465 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0])
466 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j)
467 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
468 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5]
469 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
470 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)
471 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j))
472 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
473 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2])
474 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)
475 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j))
476 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7]
477 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3])
478 * (m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
479 - m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3],
480 j)));
481 }
482
483 for(MInt j = 0; j < 3; j++) {
484 // du[j]/dx
485 c[j + 0 * 3] -=
486 F1B2 * rho * F1B2
487 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0]
488 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
489 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1], j)
490 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
491 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2]
492 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3])
493 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)
494 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2],
495 j))
496 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4]
497 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5])
498 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)
499 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4],
500 j))
501 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
502 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7])
503 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
504 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6],
505 j)));
506 // du[j]/dy
507 c[j + 1 * 3] -=
508 F1B2 * rho * F1B2
509 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2]
510 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0])
511 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2], j)
512 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
513 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3]
514 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
515 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3], j)
516 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1],
517 j))
518 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
519 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4])
520 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)
521 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4],
522 j))
523 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7]
524 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5])
525 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
526 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5],
527 j)));
528 // du[j]/dz
529 c[j + 2 * 3] -=
530 F1B2 * rho * F1B2
531 * ((m_interfaceChildren[level]->a[id].m_interpolationCoefficients[4]
532 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[0])
533 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[4], j)
534 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[0], j))
535 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[5]
536 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[1])
537 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[5], j)
538 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[1],
539 j))
540 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[6]
541 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[2])
542 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[6], j)
543 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[2],
544 j))
545 + (m_interfaceChildren[level]->a[id].m_interpolationCoefficients[7]
546 + m_interfaceChildren[level]->a[id].m_interpolationCoefficients[3])
547 * (m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[7], j)
548 - m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[3],
549 j)));
550 }
551 trace = c[0] + c[4] + c[8];
552 }
553 else {
554 cout << "Strain rate tensor not implemented for 2D yet" << endl;
555 trace = F0;
556 }
557
558 //--------------------------------------------------------------------------------------------
559
560 tmp2 = F0;
561 for(MInt d = 0; d < nDim; d++) {
562 tmp2 += (u[d] * u[d]);
563 b[2 * d] = -u[d];
564 b[2 * d + 1] = u[d];
565 }
566
567 // Calculation of new distributions
568 std::array<MFloat, nDist> eqDist;
569 eqDist = m_solver->getEqDists(rho, tmp2, u.data());
570
571 // create new incoming distributions
572 for(MInt dist = 0; dist < nDist - 1; dist++) {
573 // only the missing distributions
574 if(m_solver->a_hasNeighbor(currentId, Ld::oppositeDist(dist)) == 0) {
575 const MFloat tp = Ld::tp(Ld::distType(dist));
576
577 // put together eq-part and rescaled non-eq part; scale factor = (omega_C/omega_F) * F1B2
578
579 m_solver->a_oldDistribution(currentId, dist) = eqDist[dist] + tp * trace / omega_F;
580
581 for(MInt k = 0; k < nDim; k++) {
582 for(MInt l = 0; l < nDim; l++) {
583 m_solver->a_oldDistribution(currentId, dist) -=
584 tp * F1BCSsq * (Ld::idFld(dist, k) - 1) * (Ld::idFld(dist, l) - 1) * c[l + nDim * k] / omega_F;
585 }
586 }
587 }
588 }
589 } // end of loop over all interface cells
590 }
591 }
592 }
593}
594
601template <MInt nDim, MInt nDist, class SysEqn>
602void LbInterfaceDxQy<nDim, nDist, SysEqn>::restriction10() {
603 MFloat& tmp2 = m_static_restriction10_tmp2;
604 MFloat(&b)[2 * nDim] = m_static_restriction10_b;
605 MFloat(&c)[nDim * nDim] = m_static_restriction10_c;
606 MFloat& trace = m_static_restriction10_trace;
607
608 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
609 level++) { // level is the actual coarse resolution
610 // level+1 is the actual fine resolution
611
612 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level)
613 == 0) { // perform restriction from finest level in every second timestep,
614 // on the next coarser level every fourth timestep, and so on; starting from t = 1
615
616 for(MInt id = 0; id < m_interfaceParents[level]->size(); id++) {
617 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
618 if(currentId >= m_solver->noInternalCells()) continue;
619
620 // omega_F = 2.0 / ( 1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel()
621 // - level - 1));
622 const MFloat omega_C =
623 2.0 / (1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
624
625 MFloat rho = 0;
626 std::array<MFloat, nDim> u{};
627
628 //------------------------------------------------------------------------------------
629 // temporal extrapolation of the mean macroscopic variables and the mean stress tensor of all child cells
630 // f(t+1) = 2*f(t) - f(t-1)
631
632 for(MInt m = 0; m < IPOW2(nDim); m++) {
633 rho += 2.0 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->RHO);
634 u[0] += 2.0 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->U);
635 u[1] += 2.0 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->V);
636 IF_CONSTEXPR(nDim == 3)
637 u[2] += 2.0 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->W);
638 }
639 for(MInt m = 0; m < IPOW2(nDim); m++) {
640 rho -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->RHO);
641 u[0] -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->U);
642 u[1] -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->V);
643 IF_CONSTEXPR(nDim == 3)
644 u[2] -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->W);
645 }
646
647 // compute strain rate tensor on coarse grid: (du/dx)_coarse = 2*(du/dx)_fine
648 IF_CONSTEXPR(nDim == 3) {
649 for(MInt j = 0; j < 3; j++) {
650 // du[j]/dx
651 c[j] = 2.0 * F1B4 * rho * 2.0
652 * ((m_solver->a_variable(m_solver->c_childId(currentId, 1), j)
653 - m_solver->a_variable(m_solver->c_childId(currentId, 0), j))
654 + (m_solver->a_variable(m_solver->c_childId(currentId, 3), j)
655 - m_solver->a_variable(m_solver->c_childId(currentId, 2), j))
656 + (m_solver->a_variable(m_solver->c_childId(currentId, 5), j)
657 - m_solver->a_variable(m_solver->c_childId(currentId, 4), j))
658 + (m_solver->a_variable(m_solver->c_childId(currentId, 7), j)
659 - m_solver->a_variable(m_solver->c_childId(currentId, 6), j)));
660 // du[j]/dy
661 c[j + 3] = 2.0 * F1B4 * rho * 2.0
662 * ((m_solver->a_variable(m_solver->c_childId(currentId, 2), j)
663 - m_solver->a_variable(m_solver->c_childId(currentId, 0), j))
664 + (m_solver->a_variable(m_solver->c_childId(currentId, 3), j)
665 - m_solver->a_variable(m_solver->c_childId(currentId, 1), j))
666 + (m_solver->a_variable(m_solver->c_childId(currentId, 6), j)
667 - m_solver->a_variable(m_solver->c_childId(currentId, 4), j))
668 + (m_solver->a_variable(m_solver->c_childId(currentId, 7), j)
669 - m_solver->a_variable(m_solver->c_childId(currentId, 5), j)));
670 // du[j]/dz
671 c[j + 6] = 2.0 * F1B4 * rho * 2.0
672 * ((m_solver->a_variable(m_solver->c_childId(currentId, 4), j)
673 - m_solver->a_variable(m_solver->c_childId(currentId, 0), j))
674 + (m_solver->a_variable(m_solver->c_childId(currentId, 5), j)
675 - m_solver->a_variable(m_solver->c_childId(currentId, 1), j))
676 + (m_solver->a_variable(m_solver->c_childId(currentId, 6), j)
677 - m_solver->a_variable(m_solver->c_childId(currentId, 2), j))
678 + (m_solver->a_variable(m_solver->c_childId(currentId, 7), j)
679 - m_solver->a_variable(m_solver->c_childId(currentId, 3), j)));
680 }
681
682 for(MInt j = 0; j < 3; j++) {
683 // du[j]/dx
684 c[j] -= 1.0 * F1B4 * rho * 2.0
685 * ((m_solver->a_oldVariable(m_solver->c_childId(currentId, 1), j)
686 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 0), j))
687 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 3), j)
688 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 2), j))
689 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 5), j)
690 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 4), j))
691 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 7), j)
692 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 6), j)));
693 // du[j]/dy
694 c[j + 3] -= 1.0 * F1B4 * rho * 2.0
695 * ((m_solver->a_oldVariable(m_solver->c_childId(currentId, 2), j)
696 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 0), j))
697 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 3), j)
698 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 1), j))
699 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 6), j)
700 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 4), j))
701 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 7), j)
702 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 5), j)));
703 // du[j]/dz
704 c[j + 6] -= 1.0 * F1B4 * rho * 2.0
705 * ((m_solver->a_oldVariable(m_solver->c_childId(currentId, 4), j)
706 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 0), j))
707 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 5), j)
708 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 1), j))
709 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 6), j)
710 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 2), j))
711 + (m_solver->a_oldVariable(m_solver->c_childId(currentId, 7), j)
712 - m_solver->a_oldVariable(m_solver->c_childId(currentId, 3), j)));
713 }
714 trace = c[0] + c[4] + c[8];
715 }
716 else {
717 cout << "Strain rate tensor not implemented for 2D yet" << endl;
718 trace = F0;
719 }
720 //------------------------------------------------------------------------------------
721
722
723 // set parent values to averaged child values for subsequent prolongation
724 m_solver->a_variable(currentId, PV->RHO) = rho;
725 m_solver->a_variable(currentId, PV->U) = u[0];
726 m_solver->a_variable(currentId, PV->V) = u[1];
727 IF_CONSTEXPR(nDim == 3) m_solver->a_variable(currentId, PV->W) = u[2];
728
729 tmp2 = F0;
730 for(MInt d = 0; d < nDim; d++) {
731 tmp2 += (u[d] * u[d]);
732 b[2 * d] = -u[d];
733 b[2 * d + 1] = u[d];
734 }
735
736 std::array<MFloat, nDist> eqDist;
737 eqDist = m_solver->getEqDists(rho, tmp2, u.data());
738
739 // create new incoming distributions
740 for(MInt dist = 0; dist < nDist - 1; dist++) {
741 // only the missing distributions
742 if(!m_solver->a_hasNeighbor(currentId, Ld::oppositeDist(dist))
743 || (!m_solver->a_isInterfaceParent(m_solver->c_neighborId(currentId, Ld::oppositeDist(dist)))
744 && !m_solver->c_isLeafCell(m_solver->c_neighborId(currentId, Ld::oppositeDist(dist))))) {
745 const MFloat tp = Ld::tp(Ld::distType(dist));
746
747 // put together eq-part and rescaled non-eq part; scale factor = (omega_C/omega_F) * F1B2
748
749 m_solver->a_oldDistribution(currentId, dist) = eqDist[dist] + tp * trace / omega_C;
750
751 for(MInt k = 0; k < nDim; k++) {
752 for(MInt l = 0; l < nDim; l++) {
753 m_solver->a_oldDistribution(currentId, dist) -=
754 tp * F1BCSsq * (Ld::idFld(dist, k) - 1) * (Ld::idFld(dist, l) - 1) * c[l + nDim * k] / omega_C;
755 }
756 }
757 }
758 }
759 }
760 }
761 }
762}
763
770template <MInt nDim, MInt nDist, class SysEqn>
771template <MBool compressible>
772void LbInterfaceDxQy<nDim, nDist, SysEqn>::prolongationDupuis_() {
773#ifdef WAR_NVHPC_PSTL
774 const MInt noInternalCells = m_solver->noInternalCells();
775#endif
776
777 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
778 level++) { // minLevel+level is the actual coarse resolution
779 // minLevel+level+1 is the actual fine resolution
780
781 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1)
782 == 0) { // perform prolongation WITH temporal interpolation in every second timestep on finest level,
783 // in every fourth timestep on the next coarser level, and so on; starting from t = 2
784
785 //-----------------------------
786 // if in sync with parent:
787 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level)
788 == 0) { // perform prolongation WITHOUT temporal interpolation in every second timestep on finest level,
789 // in every fourth timestep on the next coarser level, and so on; starting from t = 1
790
791 maia::parallelFor<true>(0, m_interfaceChildren[level]->size(), [=](MInt id) {
792 const MInt currentId = m_interfaceChildren[level]->a[id].m_cellId;
793#ifdef WAR_NVHPC_PSTL
794 if(currentId >= noInternalCells) return;
795#else
796 if(currentId >= m_solver->noInternalCells()) return;
797#endif
798
799 const MFloat omega_F =
800 2.0
801 / (1.0
802 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
803 const MFloat omega_C =
804 2.0
805 / (1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
806
807 // Spatial interpolation and temporal extrapolation of macroscopic variables
808 // f(t+1) = 2*f(t) - f(t-1)
809 MFloat rho = 0.0;
810 std::array<MFloat, nDim> u{};
811
812 for(MInt m = 0; m < IPOW2(nDim); m++) {
813 rho += 2.0 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
814 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
815 for(MInt d = 0; d < nDim; d++) {
816 u[d] += 2.0 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
817 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], d);
818 }
819 }
820
821 for(MInt m = 0; m < IPOW2(nDim); m++) {
822 rho -= m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
823 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
824 for(MInt d = 0; d < nDim; d++) {
825 u[d] -= m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
826 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], d);
827 }
828 }
829
830 // Calculate all equilibrium distributions
831 std::array<MFloat, nDist> eqDist;
832 eqDist = m_solver->getEqDists(rho, u.data());
833
834 // interpolate and transform
835 for(MInt dist = 0; dist < nDist - 1; dist++) {
836 // only the missing distributions
837#ifdef WAR_NVHPC_PSTL
838 MInt oppositeDist = m_solver->m_oppositeDist[dist];
839#else
840 MInt oppositeDist = Ld::oppositeDist(dist);
841#endif
842 if(m_solver->a_hasNeighbor(currentId, dist) == 0) {
843 MFloat tmpDist1 = 0.0;
844
845 for(MInt m = 0; m < IPOW2(nDim); m++) {
846 tmpDist1 += m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
847 * m_solver->a_oldDistribution(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m],
848 oppositeDist);
849 }
850
851 m_solver->a_oldDistribution(currentId, oppositeDist) =
852 (tmpDist1 - eqDist[oppositeDist]) * (omega_C / omega_F) * F1B2 + eqDist[oppositeDist];
853 }
854 }
855 });
856 }
857
858 //---------------------------
859 // if not in sync with parent: apply temporal interpolation
860 else {
861 maia::parallelFor<true>(0, m_interfaceChildren[level]->size(), [=](MInt id) {
862 const MInt currentId = m_interfaceChildren[level]->a[id].m_cellId;
863#ifdef WAR_NVHPC_PSTL
864 if(currentId >= noInternalCells) return;
865#else
866 if(currentId >= m_solver->noInternalCells()) return;
867#endif
868
869 const MFloat omega_F =
870 2.0
871 / (1.0
872 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
873 const MFloat omega_C =
874 2.0
875 / (1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
876
877 // spatial and temporal interpolation of variables
878 // f(t+1/2) = 1.5*f(t) - 0.5*f(t-1)
879 MFloat rho = 0;
880 std::array<MFloat, nDim> u{};
881
882 for(MInt m = 0; m < IPOW2(nDim); m++) {
883 rho += 1.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
884 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
885 for(MInt d = 0; d < nDim; d++) {
886 u[d] += 1.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
887 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], d);
888 }
889 }
890
891 for(MInt m = 0; m < IPOW2(nDim); m++) {
892 rho -= 0.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
893 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
894 for(MInt d = 0; d < nDim; d++) {
895 u[d] -= 0.5 * m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
896 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], d);
897 }
898 }
899
900 // Calculation of new distributions
901 std::array<MFloat, nDist> eqDist;
902 eqDist = m_solver->getEqDists(rho, u.data());
903 std::array<MFloat, nDist> forcing{};
904 getCellForcing(forcing, currentId);
905
906 // Transform and interpolate
907 for(MInt dist = 0; dist < nDist - 1; dist++) {
908 // only for missing distributions
909#ifdef WAR_NVHPC_PSTL
910 MInt oppositeDist = m_solver->m_oppositeDist[dist];
911#else
912 MInt oppositeDist = Ld::oppositeDist(dist);
913#endif
914 if(m_solver->a_hasNeighbor(currentId, dist) == 0) {
915 MFloat tmpDist1 = 0.0;
916 MFloat tmpDist2 = 0.0;
917
918 for(MInt m = 0; m < IPOW2(nDim); m++) {
919 // add the weighted parts to the non-equilibrium distributions using linear temporal interpolation
920 // f(t+1/2) = 1/2*f(t) + 1/2*f(t+1)
921 const MFloat tmpDist = m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
922 * m_solver->a_oldDistribution(
923 m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], oppositeDist);
924 if(m_solver->a_hasNeighbor(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], dist)) {
925 // the propagation on the coarse grid has not taken place yet,
926 // thus the incoming distribution for the coarse cell can be interpolated as follows 0.5*(tmpDist1 +
927 // tmpDist2)
928 const MInt neighborId =
929 m_solver->c_neighborId(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], dist);
930 tmpDist1 += m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
931 * m_solver->a_distribution(neighborId, oppositeDist);
932 tmpDist2 += tmpDist;
933
934 // add forcing terms for coarse grid before transformation
935 // for tmpDist2 the forcing terms have already been updated
936 tmpDist1 += m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
937 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) * forcing[oppositeDist];
938 } else {
939 // if the interpolation neighbor has no neighbor in actual direction, no temporal interpolation can be
940 // applied this may only occur at a non-periodic boundary
941 tmpDist2 += tmpDist;
942
943 // forcing term is already included
944 tmpDist1 += tmpDist;
945 }
946 }
947
948 m_solver->a_oldDistribution(currentId, oppositeDist) =
949 (0.5 * (tmpDist1 + tmpDist2) - eqDist[oppositeDist]) * (omega_C / omega_F) * F1B2
950 + eqDist[oppositeDist];
951 }
952 }
953 });
954 } // end of else
955 }
956 }
957}
958
959template <MInt nDim, MInt nDist, class SysEqn>
960void LbInterfaceDxQy<nDim, nDist, SysEqn>::prolongationDupuis() {
961 prolongationDupuis_<false>();
962}
963
964template <MInt nDim, MInt nDist, class SysEqn>
965void LbInterfaceDxQy<nDim, nDist, SysEqn>::prolongationDupuisCompressible() {
966 prolongationDupuis_<true>();
967}
968
977template <MInt nDim, MInt nDist, class SysEqn>
978void LbInterfaceDxQy<nDim, nDist, SysEqn>::prolongationThermalDupuis() {
979 MFloat* ptr_m_interpolationCoefficients;
980 const MFloat factor = (m_solver->m_innerEnergy) ? ((nDim == 2) ? 3.0 : 18.0 / 5.0) : 6.0;
981
982 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
983 level++) { // minLevel+level is the actual coarse resolution
984 // minLevel+level+1 is the actual fine resolution
985
986 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1)
987 == 0) { // perform prolongation WITH temporal interpolation in every second timestep on finest level,
988 // in every fourth timestep on the next coarser level, and so on; starting from t = 2
989
990 //-----------------------------
991 // if in sync with parent:
992 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level)
993 == 0) { // perform prolongation WITHOUT temporal interpolation in every second timestep on finest level,
994 // in every fourth timestep on the next coarser level, and so on; starting from t = 1
995
996 for(MInt id = 0; id < m_interfaceChildren[level]->size(); id++) {
997 const MInt currentId = m_interfaceChildren[level]->a[id].m_cellId;
998 if(currentId >= m_solver->noInternalCells()) continue;
999
1000 ptr_m_interpolationCoefficients = m_interfaceChildren[level]->a[id].m_interpolationCoefficients;
1001
1002 const MFloat omega_F =
1003 2.0
1004 / (1.0
1005 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
1006 const MFloat omega_C =
1007 2.0
1008 / (1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
1009 const MFloat omegaT_F = 2.0
1010 / (1.0
1011 + factor * m_solver->a_kappa(currentId)
1012 * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
1013 const MFloat omegaT_C =
1014 2.0
1015 / (1.0
1016 + factor * m_solver->a_kappa(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
1017
1018 // Spatial interpolation and temporal extrapolation of macroscopic variables
1019 // f(t+1) = 2*f(t) - f(t-1)
1020 MFloat rho = 0.0;
1021 std::array<MFloat, nDim> u{};
1022 MFloat T = 0.0;
1023
1024 for(MInt m = 0; m < IPOW2(nDim); m++) {
1025 rho += 2.0 * ptr_m_interpolationCoefficients[m]
1026 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
1027 u[0] += 2.0 * ptr_m_interpolationCoefficients[m]
1028 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->U);
1029 u[1] += 2.0 * ptr_m_interpolationCoefficients[m]
1030 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->V);
1031 IF_CONSTEXPR(nDim == 3)
1032 u[2] += 2.0 * ptr_m_interpolationCoefficients[m]
1033 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->W);
1034 T += 2.0 * ptr_m_interpolationCoefficients[m]
1035 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->T);
1036 rho -= ptr_m_interpolationCoefficients[m]
1037 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
1038 u[0] -= ptr_m_interpolationCoefficients[m]
1039 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->U);
1040 u[1] -= ptr_m_interpolationCoefficients[m]
1041 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->V);
1042 IF_CONSTEXPR(nDim == 3)
1043 u[2] -= ptr_m_interpolationCoefficients[m]
1044 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->W);
1045 T -= ptr_m_interpolationCoefficients[m]
1046 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->T);
1047 }
1048
1049 // TODO labels:LB sysEqn: make PV primitive again (not storing rho*u)
1050 u[0] /= rho;
1051 u[1] /= rho;
1052 IF_CONSTEXPR(nDim == 3) u[2] /= rho;
1053
1054 std::array<MFloat, nDist> eqDist{};
1055 eqDist = m_solver->getEqDists(rho, u.data());
1056 std::array<MFloat, nDist> eqDistThermal{};
1057 eqDistThermal = m_solver->getEqDistsThermal(T, rho, u.data());
1058
1059 // interpolate and transform
1060 for(MInt dist = 0; dist < nDist - 1; dist++) {
1061 // only the missing distributions
1062 if(m_solver->a_hasNeighbor(currentId, dist) == 0) {
1063 MFloat tmpDist1 = 0.0;
1064 MFloat tmpDist1T = 0.0;
1065
1066 for(MInt m = 0; m < IPOW2(nDim); m++) {
1067 tmpDist1 += ptr_m_interpolationCoefficients[m]
1068 * m_solver->a_oldDistribution(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m],
1069 Ld::oppositeDist(dist));
1070 tmpDist1T += ptr_m_interpolationCoefficients[m]
1071 * m_solver->a_oldDistributionThermal(
1072 m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], Ld::oppositeDist(dist));
1073 }
1074
1075 m_solver->a_oldDistribution(currentId, Ld::oppositeDist(dist)) =
1076 (tmpDist1 - eqDist[Ld::oppositeDist(dist)]) * (omega_C / omega_F) * F1B2
1077 + eqDist[Ld::oppositeDist(dist)];
1078 m_solver->a_oldDistributionThermal(currentId, Ld::oppositeDist(dist)) =
1079 (tmpDist1T - eqDistThermal[Ld::oppositeDist(dist)]) * (omegaT_C / omegaT_F) * F1B2
1080 + eqDistThermal[Ld::oppositeDist(dist)];
1081 }
1082 }
1083 }
1084 }
1085
1086 //---------------------------
1087 // if not in sync with parent: apply temporal interpolation
1088 else {
1089 for(MInt id = 0; id < m_interfaceChildren[level]->size(); id++) {
1090 const MInt currentId = m_interfaceChildren[level]->a[id].m_cellId;
1091
1092 if(currentId >= m_solver->noInternalCells()) continue;
1093
1094 ptr_m_interpolationCoefficients = m_interfaceChildren[level]->a[id].m_interpolationCoefficients;
1095
1096 const MFloat omega_F =
1097 2.0
1098 / (1.0
1099 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
1100 const MFloat omega_C =
1101 2.0
1102 / (1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
1103 const MFloat omegaT_F = 2.0
1104 / (1.0
1105 + 6.0 * m_solver->a_kappa(currentId)
1106 * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
1107 const MFloat omegaT_C =
1108 2.0
1109 / (1.0
1110 + 6.0 * m_solver->a_kappa(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
1111
1112 // spatial and temporal interpolation of variables
1113 // f(t+1/2) = 1.5*f(t) - 0.5*f(t-1)
1114 MFloat rho = 0;
1115 std::array<MFloat, nDim> u{};
1116 MFloat T = 0.0;
1117
1118 for(MInt m = 0; m < IPOW2(nDim); m++) {
1119 rho += 1.5 * ptr_m_interpolationCoefficients[m]
1120 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
1121 u[0] += 1.5 * ptr_m_interpolationCoefficients[m]
1122 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->U);
1123 u[1] += 1.5 * ptr_m_interpolationCoefficients[m]
1124 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->V);
1125 IF_CONSTEXPR(nDim == 3)
1126 u[2] += 1.5 * ptr_m_interpolationCoefficients[m]
1127 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->W);
1128 T += 1.5 * ptr_m_interpolationCoefficients[m]
1129 * m_solver->a_variable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->T);
1130 rho -= 0.5 * ptr_m_interpolationCoefficients[m]
1131 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->RHO);
1132 u[0] -= 0.5 * ptr_m_interpolationCoefficients[m]
1133 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->U);
1134 u[1] -= 0.5 * ptr_m_interpolationCoefficients[m]
1135 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->V);
1136 IF_CONSTEXPR(nDim == 3)
1137 u[2] -= 0.5 * ptr_m_interpolationCoefficients[m]
1138 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->W);
1139 T -= 0.5 * ptr_m_interpolationCoefficients[m]
1140 * m_solver->a_oldVariable(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], PV->T);
1141 }
1142
1143 u[0] /= rho;
1144 u[1] /= rho;
1145 IF_CONSTEXPR(nDim == 3) u[2] /= rho;
1146
1147 std::array<MFloat, nDist> eqDist{};
1148 eqDist = m_solver->getEqDists(rho, u.data());
1149 std::array<MFloat, nDist> eqDistThermal{};
1150 eqDistThermal = m_solver->getEqDistsThermal(T, rho, u.data());
1151
1152 std::array<MFloat, nDist> forcing{};
1153 getCellForcing(forcing, currentId);
1154
1155 // Transform and interpolate
1156 for(MInt dist = 0; dist < nDist - 1; dist++) {
1157 // only for missing distributions
1158 if(m_solver->a_hasNeighbor(currentId, dist) == 0) {
1159 MFloat tmpDist1 = 0.0;
1160 MFloat tmpDist2 = 0.0;
1161 MFloat tmpDist1T = 0.0;
1162 MFloat tmpDist2T = 0.0;
1163
1164 for(MInt m = 0; m < IPOW2(nDim); m++) {
1165 // add the weighted parts to the non-equilibrium distributions using linear temporal interpolation
1166 // f(t+1/2) = 1/2*f(t) + 1/2*f(t+1)
1167 if(m_solver->a_hasNeighbor(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], dist)) {
1168 const MInt neighborId =
1169 m_solver->c_neighborId(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], dist);
1170
1171 // the propagation on the coarse grid has not taken place yet,
1172 // thus the incoming distribution for the coarse cell can be interpolated as follows 0.5*(tmpDist1 +
1173 // tmpDist2)
1174 tmpDist1 +=
1175 ptr_m_interpolationCoefficients[m] * m_solver->a_distribution(neighborId, Ld::oppositeDist(dist));
1176 tmpDist2 +=
1177 ptr_m_interpolationCoefficients[m]
1178 * m_solver->a_oldDistribution(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m],
1179 Ld::oppositeDist(dist));
1180 tmpDist1T += ptr_m_interpolationCoefficients[m]
1181 * m_solver->a_distributionThermal(neighborId, Ld::oppositeDist(dist));
1182 tmpDist2T +=
1183 ptr_m_interpolationCoefficients[m]
1184 * m_solver->a_oldDistributionThermal(
1185 m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], Ld::oppositeDist(dist));
1186
1187 // add forcing terms for coarse grid before transformation
1188 // for tmpDist2 the forcing terms have already been updated
1189 tmpDist1 += m_interfaceChildren[level]->a[id].m_interpolationCoefficients[m]
1190 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level)
1191 * forcing[Ld::oppositeDist(dist)];
1192 }
1193
1194 // if the interpolation neighbor has no neighbor in actual direction, no temporal interpolation can be
1195 // applied this may only occur at a non-periodic boundary
1196 else {
1197 tmpDist2 +=
1198 ptr_m_interpolationCoefficients[m]
1199 * m_solver->a_oldDistribution(m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m],
1200 Ld::oppositeDist(dist));
1201 tmpDist2T +=
1202 ptr_m_interpolationCoefficients[m]
1203 * m_solver->a_oldDistributionThermal(
1204 m_interfaceChildren[level]->a[id].m_interpolationNeighbors[m], Ld::oppositeDist(dist));
1205
1206 // set tmpDist2=tmpDist1
1207 // forcing term is already included
1208 tmpDist1 = tmpDist2;
1209 tmpDist1T = tmpDist2T;
1210 }
1211 }
1212
1213 m_solver->a_oldDistribution(currentId, Ld::oppositeDist(dist)) =
1214 (0.5 * (tmpDist1 + tmpDist2) - eqDist[Ld::oppositeDist(dist)]) * (omega_C / omega_F) * F1B2
1215 + eqDist[Ld::oppositeDist(dist)];
1216 m_solver->a_oldDistributionThermal(currentId, Ld::oppositeDist(dist)) =
1217 (0.5 * (tmpDist1T + tmpDist2T) - eqDistThermal[Ld::oppositeDist(dist)]) * (omegaT_C / omegaT_F) * F1B2
1218 + eqDistThermal[Ld::oppositeDist(dist)];
1219 }
1220 }
1221 }
1222 } // end of else
1223 }
1224 }
1225}
1226
1233template <MInt nDim, MInt nDist, class SysEqn>
1234template <MBool compressible>
1235void LbInterfaceDxQy<nDim, nDist, SysEqn>::restrictionDupuis_() {
1236 TRACE();
1237
1238#ifdef WAR_NVHPC_PSTL
1239 const MInt noInternalCells = m_solver->noInternalCells();
1240#endif
1241
1242 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
1243 level++) { // level is the actual coarse resolution
1244 // level+1 is the actual fine resolution
1245
1246 const MFloat nu_F = m_solver->m_Ma * LBCS / m_solver->m_Re * m_solver->m_referenceLength
1247 * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1);
1248 const MFloat omega_F = 2.0 / (1.0 + 6.0 * nu_F);
1249 const MFloat nu_C = m_solver->m_Ma * LBCS / m_solver->m_Re * m_solver->m_referenceLength
1250 * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level);
1251 const MFloat omega_C = 2.0 / (1.0 + 6.0 * nu_C);
1252
1253
1254 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level)
1255 == 0) { // perform restriction from finest level in every second timestep,
1256 // on the next coarser level every fourth timestep, and so on; starting from t = 1
1257
1258 maia::parallelFor<true>(0, m_interfaceParents[level]->size(), [=](MInt id) {
1259 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
1260#ifdef WAR_NVHPC_PSTL
1261 if(currentId >= noInternalCells) return;
1262#else
1263 if(currentId >= m_solver->noInternalCells()) return;
1264#endif
1265
1266 MFloat rho = 0;
1267 std::array<MFloat, nDim> u{};
1268
1269 // temporal extrapolated mean of the macroscopic variables of the child cells
1270 // f(t+1) = 2*f(t) - f(t-1)
1271 for(MInt m = 0; m < IPOW2(nDim); m++) {
1272 const MInt childCellId = m_solver->c_childId(currentId, m);
1273 rho += FFPOW2(nDim)
1274 * (2.0 * m_solver->a_variable(childCellId, PV->RHO) - m_solver->a_oldVariable(childCellId, PV->RHO));
1275 for(MInt n = 0; n < nDim; n++) {
1276 u[n] +=
1277 FFPOW2(nDim) * (2.0 * m_solver->a_variable(childCellId, n) - m_solver->a_oldVariable(childCellId, n));
1278 }
1279 }
1280
1281 std::array<MFloat, nDist> eqDist;
1282 eqDist = m_solver->getEqDists(rho, u.data());
1283
1284 // Interpolate missing distributions
1285 for(MInt dist = 0; dist < (nDist - 1); dist++) {
1286#ifdef WAR_NVHPC_PSTL
1287 MInt oppositeDist = m_solver->m_oppositeDist[dist];
1288#else
1289 MInt oppositeDist = Ld::oppositeDist(dist);
1290#endif
1291 if(!m_solver->a_isInterfaceParent(m_solver->c_neighborId(currentId, oppositeDist))
1292 && !m_solver->c_isLeafCell(m_solver->c_neighborId(currentId, oppositeDist))) {
1293 // Calculate mean distribution from all children
1294 MFloat tmpDist = 0.0;
1295 for(MInt m = 0; m < IPOW2(nDim); m++) {
1296 tmpDist += FFPOW2(nDim) * m_solver->a_oldDistribution(m_solver->c_childId(currentId, m), dist);
1297 }
1298
1299 // transform distribution
1300 m_solver->a_oldDistribution(currentId, dist) =
1301 (tmpDist - eqDist[dist]) * (omega_F / omega_C) * 2.0 + eqDist[dist];
1302 }
1303 }
1304 });
1305 }
1306 }
1307}
1308
1309template <MInt nDim, MInt nDist, class SysEqn>
1310void LbInterfaceDxQy<nDim, nDist, SysEqn>::restrictionDupuis() {
1311 restrictionDupuis_<false>();
1312}
1313
1314template <MInt nDim, MInt nDist, class SysEqn>
1315void LbInterfaceDxQy<nDim, nDist, SysEqn>::restrictionDupuisCompressible() {
1316 restrictionDupuis_<true>();
1317}
1318
1327template <MInt nDim, MInt nDist, class SysEqn>
1328void LbInterfaceDxQy<nDim, nDist, SysEqn>::restrictionThermalDupuis() {
1329 const MFloat factor = (m_solver->m_innerEnergy) ? ((nDim == 2) ? 3.0 : 18.0 / 5.0) : 6.0;
1330 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
1331 level++) { // level is the actual coarse resolution
1332 // level+1 is the actual fine resolution
1333
1334 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level)
1335 == 0) { // perform restriction from finest level in every second timestep,
1336 // on the next coarser level every fourth timestep, and so on; starting from t = 1
1337
1338 for(MInt id = 0; id < m_interfaceParents[level]->size(); id++) {
1339 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
1340 if(currentId >= m_solver->noInternalCells()) continue;
1341
1342 const MFloat omega_F =
1343 2.0
1344 / (1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
1345 const MFloat omega_C =
1346 2.0 / (1.0 + 6.0 * m_solver->a_nu(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
1347 const MFloat omegaT_F = 2.0
1348 / (1.0
1349 + factor * m_solver->a_kappa(currentId)
1350 * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1));
1351 const MFloat omegaT_C =
1352 2.0
1353 / (1.0
1354 + factor * m_solver->a_kappa(currentId) * FFPOW2(m_solver->maxLevel() - m_solver->minLevel() - level));
1355
1356 MFloat rho = F0;
1357 std::array<MFloat, nDim> u{};
1358 MFloat T = F0;
1359
1360 // temporal extrapolated mean of the macroscopic variables of the child cells
1361 // f(t+1) = 2*f(t) - f(t-1)
1362 for(MInt m = 0; m < IPOW2(nDim); m++) {
1363 rho += F2 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->RHO);
1364 u[0] += F2 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->U);
1365 u[1] += F2 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->V);
1366 IF_CONSTEXPR(nDim == 3)
1367 u[2] += F2 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->W);
1368 T += F2 * FFPOW2(nDim) * m_solver->a_variable(m_solver->c_childId(currentId, m), PV->T);
1369 rho -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->RHO);
1370 u[0] -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->U);
1371 u[1] -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->V);
1372 IF_CONSTEXPR(nDim == 3)
1373 u[2] -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->W);
1374 T -= FFPOW2(nDim) * m_solver->a_oldVariable(m_solver->c_childId(currentId, m), PV->T);
1375 }
1376
1377 for(MInt d = 0; d < nDim; d++) {
1378 u[d] /= rho;
1379 }
1380
1381 std::array<MFloat, nDist> eqDist{};
1382 eqDist = m_solver->getEqDists(rho, u.data());
1383 std::array<MFloat, nDist> eqDistThermal{};
1384 eqDistThermal = m_solver->getEqDistsThermal(T, rho, u.data());
1385
1386 // Interpolate missing distributions
1387 for(MInt dist = 0; dist < (nDist - 1); dist++) {
1388 if(!m_solver->a_isInterfaceParent(m_solver->c_neighborId(currentId, Ld::oppositeDist(dist)))
1389 && !m_solver->c_isLeafCell(m_solver->c_neighborId(currentId, Ld::oppositeDist(dist)))) {
1390 // Calculate mean distribution from all children
1391 MFloat tmpDist = 0.0;
1392 MFloat tmpDistT = 0.0;
1393 for(MInt m = 0; m < IPOW2(nDim); m++) {
1394 tmpDist += FFPOW2(nDim) * m_solver->a_oldDistribution(m_solver->c_childId(currentId, m), dist);
1395 tmpDistT += FFPOW2(nDim) * m_solver->a_oldDistributionThermal(m_solver->c_childId(currentId, m), dist);
1396 }
1397
1398 // transform distribution
1399 m_solver->a_oldDistribution(currentId, dist) =
1400 (tmpDist - eqDist[dist]) * (omega_F / omega_C) * 2.0 + eqDist[dist];
1401 m_solver->a_oldDistributionThermal(currentId, dist) =
1402 (tmpDistT - eqDistThermal[dist]) * (omegaT_F / omegaT_C) * 2.0 + eqDistThermal[dist];
1403 }
1404 }
1405 }
1406 }
1407 }
1408}
1409
1415template <MInt nDim, MInt nDist, class SysEqn>
1416void LbInterfaceDxQy<nDim, nDist, SysEqn>::prolongationRohde() {
1417 std::array<MFloat, nDist> forcing{};
1418 getCellForcing(forcing);
1419
1420 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
1421 level++) { // minLevel+level is the actual coarse resolution
1422 // minLevel+level+1 is the actual fine resolution
1423
1424 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1) == 0) {
1425 //-----------------------
1426 // intermittant prolongation
1427 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) == 0) {
1428 for(MInt id = 0; id < m_interfaceParents[level]->size(); id++) {
1429 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
1430 if(currentId >= m_solver->noInternalCells()) continue;
1431 for(MInt m = 0; m < IPOW2(nDim); m++) {
1432 // compensate collision on fine grid
1433 if(m_cellDependentForcing) getCellForcing(forcing, m_solver->c_childId(currentId, m));
1434 for(MInt dist = 0; dist < nDist; dist++) {
1435 // remove forcing term
1436 m_solver->a_oldDistribution(m_solver->c_childId(currentId, m), dist) -=
1437 FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1) * forcing[dist];
1438
1439 // reset outgoing distributions
1440 m_solver->a_distribution(m_solver->c_childId(currentId, m), dist) =
1441 m_solver->a_oldDistribution(m_solver->c_childId(currentId, m), dist);
1442 }
1443 // copy variables from parent
1444 for(MInt var = 0; var < (nDim + 1); var++) {
1445 m_solver->a_variable(m_solver->c_childId(currentId, m), var) = m_solver->a_variable(currentId, var);
1446 }
1447 }
1448 }
1449
1450 }
1451 //--------------------------
1452 // regular prolongation
1453 else {
1454 for(MInt id = 0; id < m_interfaceParents[level]->size(); id++) {
1455 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
1456 if(currentId >= m_solver->noInternalCells()) continue;
1457
1458 // overwrite variables and distributions
1459 for(MInt m = 0; m < IPOW2(nDim); m++) {
1460 // hand over distributions
1461 if(m_cellDependentForcing) getCellForcing(forcing, m_solver->c_childId(currentId, m));
1462 for(MInt dist = 0; dist < nDist; dist++) {
1463 m_solver->a_distribution(m_solver->c_childId(currentId, m), dist) =
1464 m_solver->a_distribution(currentId, dist);
1465
1466 // add half of forcing term for coarse grid
1467 m_solver->a_distribution(m_solver->c_childId(currentId, m), dist) +=
1468 F1B2 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) * forcing[dist];
1469 // remove half of forcing term for fine grid
1470 m_solver->a_distribution(m_solver->c_childId(currentId, m), dist) -=
1471 F1B2 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1) * forcing[dist];
1472 }
1473
1474 // copy variables from parent
1475 for(MInt var = 0; var < (nDim + 1); var++) {
1476 m_solver->a_variable(m_solver->c_childId(currentId, m), var) = m_solver->a_variable(currentId, var);
1477 }
1478 }
1479 }
1480 }
1481 }
1482 }
1483}
1484
1493template <MInt nDim, MInt nDist, class SysEqn>
1494void LbInterfaceDxQy<nDim, nDist, SysEqn>::prolongationThermalRohde() {
1495 std::array<MFloat, nDist> forcing{};
1496 getCellForcing(forcing);
1497
1498 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
1499 level++) { // minLevel+level is the actual coarse resolution
1500 // minLevel+level+1 is the actual fine resolution
1501
1502 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1) == 0) {
1503 //-----------------------
1504 // intermittant prolongation
1505 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) == 0) {
1506 for(MInt id = 0; id < m_interfaceParents[level]->size(); id++) {
1507 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
1508 if(currentId >= m_solver->noInternalCells()) continue;
1509
1510 for(MInt m = 0; m < IPOW2(nDim); m++) {
1511 // compensate collision on fine grid
1512 if(m_cellDependentForcing) getCellForcing(forcing, m_solver->c_childId(currentId, m));
1513 for(MInt dist = 0; dist < nDist; dist++) {
1514 // remove forcing term
1515 m_solver->a_oldDistribution(m_solver->c_childId(currentId, m), dist) -=
1516 FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1) * forcing[dist];
1517
1518 // reset outgoing distributions
1519 m_solver->a_distribution(m_solver->c_childId(currentId, m), dist) =
1520 m_solver->a_oldDistribution(m_solver->c_childId(currentId, m), dist);
1521 m_solver->a_distributionThermal(m_solver->c_childId(currentId, m), dist) =
1522 m_solver->a_oldDistributionThermal(m_solver->c_childId(currentId, m), dist);
1523 }
1524 // copy variables from parent
1525 for(MInt var = 0; var < (nDim + 2); var++) {
1526 m_solver->a_variable(m_solver->c_childId(currentId, m), var) = m_solver->a_variable(currentId, var);
1527 }
1528 }
1529 }
1530
1531 }
1532 //--------------------------
1533 // regular prolongation
1534 else {
1535 for(MInt id = 0; id < m_interfaceParents[level]->size(); id++) {
1536 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
1537 if(currentId >= m_solver->noInternalCells()) continue;
1538
1539 // overwrite variables and distributions
1540 for(MInt m = 0; m < IPOW2(nDim); m++) {
1541 // hand over distributions
1542 if(m_cellDependentForcing) getCellForcing(forcing, m_solver->c_childId(currentId, m));
1543 for(MInt dist = 0; dist < nDist; dist++) {
1544 m_solver->a_distribution(m_solver->c_childId(currentId, m), dist) =
1545 m_solver->a_distribution(currentId, dist);
1546 m_solver->a_distributionThermal(m_solver->c_childId(currentId, m), dist) =
1547 m_solver->a_distributionThermal(currentId, dist);
1548
1549 // add half of forcing term for coarse grid
1550 m_solver->a_distribution(m_solver->c_childId(currentId, m), dist) +=
1551 F1B2 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) * forcing[dist];
1552 // remove half of forcing term for fine grid
1553 m_solver->a_distribution(m_solver->c_childId(currentId, m), dist) -=
1554 F1B2 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1) * forcing[dist];
1555 }
1556
1557 // copy variables from parent
1558 for(MInt var = 0; var < (nDim + 2); var++) {
1559 m_solver->a_variable(m_solver->c_childId(currentId, m), var) = m_solver->a_variable(currentId, var);
1560 }
1561 }
1562 }
1563 }
1564 }
1565 }
1566}
1567
1568
1575template <MInt nDim, MInt nDist, class SysEqn>
1576void LbInterfaceDxQy<nDim, nDist, SysEqn>::restrictionRohde() {
1577 std::array<MFloat, nDist> forcing{};
1578 getCellForcing(forcing);
1579
1580 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
1581 level++) { // minLevel+level is the actual coarse resolution
1582 // minLevel+level+1 is the actual fine resolution
1583
1584 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) == 0) {
1585 for(MInt id = 0; id < m_interfaceParents[level]->size(); id++) {
1586 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
1587 if(currentId >= m_solver->noInternalCells()) continue;
1588 if(m_cellDependentForcing) getCellForcing(forcing, currentId);
1589
1590 for(MInt dist = 0; dist < (nDist - 1); dist++) {
1591 if(!m_solver->a_isInterfaceParent(m_solver->c_neighborId(currentId, Ld::oppositeDist(dist)))
1592 && !m_solver->c_isLeafCell(m_solver->c_neighborId(currentId, Ld::oppositeDist(dist)))) {
1593 m_solver->a_oldDistribution(currentId, dist) = F0;
1594
1595 for(MInt m = 0; m < IPOW2(nDim); m++) {
1596 m_solver->a_oldDistribution(currentId, dist) +=
1597 FFPOW2(nDim) * m_solver->a_oldDistribution(m_solver->c_childId(currentId, m), dist);
1598 }
1599
1600 // add half of forcing term for fine grid
1601 m_solver->a_oldDistribution(currentId, dist) +=
1602 F1B2 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1) * forcing[dist];
1603 // compensate half of forcing term for coarse grid
1604 m_solver->a_oldDistribution(currentId, dist) -=
1605 F1B2 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) * forcing[dist];
1606 }
1607 }
1608 }
1609 }
1610 }
1611}
1612
1621template <MInt nDim, MInt nDist, class SysEqn>
1622void LbInterfaceDxQy<nDim, nDist, SysEqn>::restrictionThermalRohde() {
1623 std::array<MFloat, nDist> forcing{};
1624 getCellForcing(forcing);
1625
1626 for(MInt level = 0; level < (m_solver->maxLevel() - m_solver->minLevel());
1627 level++) { // minLevel+level is the actual coarse resolution
1628 // minLevel+level+1 is the actual fine resolution
1629
1630 if((globalTimeStep) % IPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) == 0) {
1631 for(MInt id = 0; id < m_interfaceParents[level]->size(); id++) {
1632 const MInt currentId = m_interfaceParents[level]->a[id].m_cellId;
1633 if(currentId >= m_solver->noInternalCells()) continue;
1634 if(m_cellDependentForcing) getCellForcing(forcing, currentId);
1635
1636 for(MInt dist = 0; dist < (nDist - 1); dist++) {
1637 if(!m_solver->a_isInterfaceParent(m_solver->c_neighborId(currentId, Ld::oppositeDist(dist)))
1638 && !m_solver->c_isLeafCell(m_solver->c_neighborId(currentId, Ld::oppositeDist(dist)))) {
1639 m_solver->a_oldDistribution(currentId, dist) = F0;
1640 m_solver->a_oldDistributionThermal(currentId, dist) = F0;
1641
1642 for(MInt m = 0; m < IPOW2(nDim); m++) {
1643 m_solver->a_oldDistribution(currentId, dist) +=
1644 FFPOW2(nDim) * m_solver->a_oldDistribution(m_solver->c_childId(currentId, m), dist);
1645 m_solver->a_oldDistributionThermal(currentId, dist) +=
1646 FFPOW2(nDim) * m_solver->a_oldDistributionThermal(m_solver->c_childId(currentId, m), dist);
1647 }
1648
1649 // add half of forcing term for fine grid
1650 m_solver->a_oldDistribution(currentId, dist) +=
1651 F1B2 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level - 1) * forcing[dist];
1652 // compensate half of forcing term for coarse grid
1653 m_solver->a_oldDistribution(currentId, dist) -=
1654 F1B2 * FPOW2(m_solver->maxLevel() - m_solver->minLevel() - level) * forcing[dist];
1655 }
1656 }
1657 }
1658 }
1659 }
1660}
1661
1662//------------------ ADAPTATION -----------------------//
1663
1669template <MInt nDim, MInt nDist, class SysEqn>
1670void LbInterfaceDxQy<nDim, nDist, SysEqn>::setAdaptationFunctions() {
1671 TRACE();
1672
1673 switch(string2enum(m_adaptationInitMethod)) {
1674 // Default property value
1675 case INIT_DUPUIS_FILIPPOVA: {
1676 fRefineCell = &LbInterfaceDxQy::refineCellDupuis;
1677 fRemoveChildren = &LbInterfaceDxQy::removeChildsDupuisFilippova;
1678 break;
1679 }
1680 case INIT_COPYPASTE: {
1681 fRefineCell = &LbInterfaceDxQy::refineCellCopyPaste;
1682 fRemoveChildren = &LbInterfaceDxQy::removeChildsCopyPaste;
1683 break;
1684 }
1685 default: {
1686 stringstream errorMessage;
1687 errorMessage << " LbInterface::setAdaptationFunctions: Specified interface condition: " << m_adaptationInitMethod
1688 << " does not exist. Exiting!";
1689 mTerm(1, AT_, errorMessage.str());
1690 }
1691 }
1692}
1693
1694template <MInt nDim, MInt nDist, class SysEqn>
1695void LbInterfaceDxQy<nDim, nDist, SysEqn>::refineCell(const MInt parentId, const MInt* childIds) {
1696 (this->*fRefineCell)(parentId, childIds);
1697}
1698
1699template <MInt nDim, MInt nDist, class SysEqn>
1700void LbInterfaceDxQy<nDim, nDist, SysEqn>::removeChildren(const MInt parentId) {
1701 (this->*fRemoveChildren)(parentId);
1702}
1703
1713template <MInt nDim, MInt nDist, class SysEqn>
1714void LbInterfaceDxQy<nDim, nDist, SysEqn>::refineCellDupuis(const MInt parentId, const MInt* childIds) {
1715 // Halo cells are initialized by exchange
1716 // ToDo: This check is unnecessary by now
1717 // because this function should not be called for Halos anymore.
1718 if(!m_solver->a_isHalo(parentId)) {
1719 // Find parent neighbors
1720 // ToDo: Check if c_neighbor of solver can be used. Should be possible now because
1721 // this function is calles at the end of reinitAfterAdaptation and c_neighbor is
1722 // updated before.
1723 MInt parentNeighbors[IPOW3[nDim] - 1];
1724 for(MInt n = 0; n < IPOW3[nDim] - 1; n++) {
1725 parentNeighbors[n] = m_solver->c_neighborId(parentId, n);
1726 }
1727
1728 const MInt parentLevel = m_solver->c_level(parentId);
1729 const MInt dParentLevel = m_solver->maxLevel() - parentLevel;
1730 const MFloat omega_F = 2.0 / (1.0 + 6.0 * m_solver->a_nu(parentId) * FFPOW2(dParentLevel - 1));
1731 const MFloat omega_C = 2.0 / (1.0 + 6.0 * m_solver->a_nu(parentId) * FFPOW2(dParentLevel));
1732 std::array<MFloat, nDist> forcing{};
1733 getCellForcing(forcing);
1734
1735 // Parent and child level have done propagation in the last timestep and will do collision after this...
1736 if((globalTimeStep - 1) % IPOW2(dParentLevel) == 0) {
1737 // Loop over all children that are created by refining the cell
1738 for(MInt child = 0; child < IPOW2(nDim); child++) {
1739 if(childIds[child] < 0) continue;
1740 const MInt childId = childIds[child];
1741
1742 m_solver->a_nu(childId) = m_solver->a_nu(parentId);
1743
1744 // Get position of child in parent
1745 MInt position = 0;
1746 for(MInt dim = 0; dim < nDim; dim++) {
1747 if(m_solver->a_coordinate(childId, dim) < m_solver->a_coordinate(parentId, dim)) {
1748 position += 0;
1749 } else {
1750 position += IPOW2(dim);
1751 }
1752 }
1753
1754 // Set interpolation neighbors
1755 // All interpolation neighbors are on parent level
1756 MInt interpolationNeighbors[IPOW2(nDim)];
1757 for(MInt k = 0; k < IPOW2(nDim); k++) {
1758 const MInt dir = Ld::intNghbrArray(position, k);
1759 if(dir < IPOW3[nDim] - 1) {
1760 interpolationNeighbors[k] = parentNeighbors[dir];
1761 } else {
1762 interpolationNeighbors[k] = parentId;
1763 }
1764 }
1765 // Set interpolation coefficients
1766 MFloat interpolationCoefficients[IPOW2(nDim)];
1767 for(MInt k = 0; k < IPOW2(nDim); k++) {
1768 interpolationCoefficients[k] = Ld::linearInterpolationCoefficients(position, k);
1769 }
1770
1771 MFloat rho = 0;
1772 std::array<MFloat, nDim> u{};
1773 MFloat rhoOld = 0;
1774 std::array<MFloat, nDim> uOld{};
1775
1776 // Spatial interpolation
1777 for(MInt m = 0; m < IPOW2(nDim); m++) {
1778 const MInt interpNeighbor = interpolationNeighbors[m];
1779 const MFloat interpCoeff = interpolationCoefficients[m];
1780
1781 // Macroscopic vars at t
1782 rho += interpCoeff * m_solver->a_variable(interpNeighbor, PV->RHO);
1783 for(MInt n = 0; n < nDim; n++) {
1784 u[n] += interpCoeff * m_solver->a_variable(interpNeighbor, PV->VV[n]);
1785 }
1786
1787 // Macroscopic vars at t-1
1788 rhoOld += interpCoeff * m_solver->a_variable(interpNeighbor, PV->RHO);
1789 for(MInt n = 0; n < nDim; n++) {
1790 uOld[n] += interpCoeff * m_solver->a_variable(interpNeighbor, PV->VV[n]);
1791 }
1792 }
1793
1794 // Temporal extrapolation: f(t+1) = 2*f(t) - f(t-1)
1795 // Needed for initialization of oldDistribution
1796 const MFloat rho_tp1 = 2 * rho - rhoOld;
1797 std::array<MFloat, nDim> u_tp1{};
1798 for(MInt d = 0; d < nDim; d++) {
1799 u_tp1[d] = 2 * u[d] - uOld[d];
1800 }
1801
1802 // Temporal extrapolation: f(t+1/2) = 1.5*f(t) - 0.5*f(t-1)
1803 // Needed for initialization of macroscopic variables
1804 MFloat rho_tp12 = 1.5 * rho - 0.5 * rhoOld;
1805 std::array<MFloat, nDim> u_tp12{};
1806 for(MInt d = 0; d < nDim; d++) {
1807 u_tp12[d] = 1.5 * u[d] - 0.5 * uOld[d];
1808 }
1809
1810 // Setting macroscopic variables in child cell
1811 m_solver->a_variable(childId, PV->RHO) = rho_tp12;
1812 for(MInt n = 0; n < nDim; n++) {
1813 m_solver->a_variable(childId, PV->VV[n]) = u_tp12[n];
1814 }
1815 m_solver->a_oldVariable(childId, PV->RHO) = rho;
1816 for(MInt n = 0; n < nDim; n++) {
1817 m_solver->a_oldVariable(childId, PV->VV[n]) = u[n];
1818 }
1819
1820 // Calculate all equilibrium distribution for tp1
1821 std::array<MFloat, nDist> eqDist_tp1{};
1822 eqDist_tp1 = m_solver->getEqDists(rho_tp1, u_tp1.data());
1823
1824 // Transform and interpolate all distributions
1825 for(MInt dist = 0; dist < nDist; dist++) {
1826 MFloat tmpDist = 0.0;
1827 for(MInt m = 0; m < IPOW2(nDim); m++) {
1828 tmpDist += interpolationCoefficients[m] * m_solver->a_oldDistribution(interpolationNeighbors[m], dist);
1829 }
1830
1831 // Distribution before next collision
1832 m_solver->a_oldDistribution(childId, dist) =
1833 (tmpDist - eqDist_tp1[dist]) * (omega_C / omega_F) * F1B2 + eqDist_tp1[dist];
1834
1835 // All cells need to perform collision step next to initialize this variable
1836 m_solver->a_distribution(childId, dist) = m_solver->a_oldDistribution(childId, dist);
1837 }
1838 }
1839 } else {
1840 // Loop over all children that are created by refining the cell
1841 for(MInt child = 0; child < IPOW2(nDim); child++) {
1842 if(childIds[child] < 0) continue;
1843 const MInt childId = childIds[child];
1844
1845 m_solver->a_nu(childId) = m_solver->a_nu(parentId);
1846
1847 // Get position of child in parent
1848 MInt position = 0;
1849 for(MInt dim = 0; dim < nDim; dim++) {
1850 if(m_solver->a_coordinate(childId, dim) < m_solver->a_coordinate(parentId, dim)) {
1851 position += 0;
1852 } else {
1853 position += IPOW2(dim);
1854 }
1855 }
1856
1857 // Set interpolation neighbors
1858 // All interpolation neighbors are on parent level
1859 MInt interpolationNeighbors[IPOW2(nDim)];
1860 for(MInt k = 0; k < IPOW2(nDim); k++) {
1861 const MInt dir = Ld::intNghbrArray(position, k);
1862 if(dir < IPOW3[nDim] - 1) {
1863 interpolationNeighbors[k] = parentNeighbors[dir];
1864 } else {
1865 interpolationNeighbors[k] = parentId;
1866 }
1867 }
1868 // Set interpolation coefficients
1869 MFloat interpolationCoefficients[IPOW2(nDim)];
1870 for(MInt k = 0; k < IPOW2(nDim); k++) {
1871 interpolationCoefficients[k] = Ld::linearInterpolationCoefficients(position, k);
1872 }
1873
1874 MFloat rho = 0;
1875 std::array<MFloat, nDim> u{};
1876 MFloat rhoOld = 0;
1877 std::array<MFloat, nDim> uOld{};
1878
1879 // Spatial interpolation
1880 for(MInt m = 0; m < IPOW2(nDim); m++) {
1881 const MInt interpNeighbor = interpolationNeighbors[m];
1882 const MFloat interpCoeff = interpolationCoefficients[m];
1883
1884 // Macroscopic vars at t
1885 rho += interpCoeff * m_solver->a_variable(interpNeighbor, PV->RHO);
1886 for(MInt n = 0; n < nDim; n++) {
1887 u[n] += interpCoeff * m_solver->a_variable(interpNeighbor, PV->VV[n]);
1888 }
1889
1890 // Macroscopic vars at t-1
1891 rhoOld += interpCoeff * m_solver->a_variable(interpNeighbor, PV->RHO);
1892 for(MInt n = 0; n < nDim; n++) {
1893 uOld[n] += interpCoeff * m_solver->a_variable(interpNeighbor, PV->VV[n]);
1894 }
1895 }
1896
1897 // Temporal extrapolation: f(t+1/2) = 1.5*f(t) - 0.5*f(t-1)
1898 // Needed for initialization of oldDistribution
1899 const MFloat rho_tp1 = 1.5 * rho - 0.5 * rhoOld;
1900 std::array<MFloat, nDim> u_tp1{};
1901 for(MInt d = 0; d < nDim; d++) {
1902 u_tp1[d] = 1.5 * u[d] - 0.5 * uOld[d];
1903 }
1904
1905 // Temporal extrapolation: f(t-1/2) = 0.5*f(t) + 0.5*f(t-1)
1906 // Needed for initialization of macroscopic variables
1907 MFloat rho_tp12 = 0.5 * rho + 0.5 * rhoOld;
1908 std::array<MFloat, nDim> u_tp12{};
1909 for(MInt d = 0; d < nDim; d++) {
1910 u_tp12[d] = 0.5 * u[d] + 0.5 * uOld[d];
1911 }
1912
1913 // Setting macroscopic variables in child cell
1914 m_solver->a_variable(childId, PV->RHO) = rho;
1915 for(MInt n = 0; n < nDim; n++) {
1916 m_solver->a_variable(childId, PV->VV[n]) = u[n];
1917 }
1918 m_solver->a_oldVariable(childId, PV->RHO) = rho_tp12;
1919 for(MInt n = 0; n < nDim; n++) {
1920 m_solver->a_oldVariable(childId, PV->VV[n]) = u_tp12[n];
1921 }
1922
1923 // Calculate all equilibrium distribution for tp1
1924 std::array<MFloat, nDist> eqDist_tp1{};
1925 eqDist_tp1 = m_solver->getEqDists(rho_tp1, u_tp1.data());
1926
1927 // Transform and interpolate all distributions
1928 MFloat sumOfOldDist = F0;
1929 for(MInt dist = 0; dist < nDist - 1; dist++) {
1930 MBool propagateNeighborValue = false;
1931
1932 if(m_solver->a_hasNeighbor(childId, Ld::oppositeDist(dist))) {
1933 if(!m_solver->a_hasProperty(childId, LbCell::WasNewlyCreated)) {
1934 propagateNeighborValue = true;
1935 }
1936 }
1937
1938 if(propagateNeighborValue) {
1939 m_solver->a_oldDistribution(childId, dist) =
1940 m_solver->a_distribution(m_solver->c_neighborId(childId, Ld::oppositeDist(dist)), dist)
1941 + FPOW2(m_solver->maxLevel() - m_solver->a_level(childId)) * forcing[dist];
1942 } else {
1943 MFloat tmpDist1 = 0.0;
1944 MFloat tmpDist2 = 0.0;
1945 for(MInt m = 0; m < IPOW2(nDim); m++) {
1946 MFloat tmpDist =
1947 interpolationCoefficients[m] * m_solver->a_oldDistribution(interpolationNeighbors[m], dist);
1948 if(m_solver->a_hasNeighbor(interpolationNeighbors[m], Ld::oppositeDist(dist))) {
1949 const MInt neighborId = m_solver->c_neighborId(interpolationNeighbors[m], Ld::oppositeDist(dist));
1950 tmpDist1 += interpolationCoefficients[m] * m_solver->a_distribution(neighborId, dist);
1951 tmpDist1 += interpolationCoefficients[m] * FPOW2(m_solver->maxLevel() - m_solver->a_level(childId) + 1)
1952 * forcing[dist];
1953 tmpDist2 += tmpDist;
1954 } else {
1955 // if the interpolation neighbor has no neighbor in actual direction, no temporal interpolation can be
1956 // applied this may only occur at a non-periodic boundary
1957 tmpDist2 += tmpDist;
1958 tmpDist1 += tmpDist;
1959 }
1960 }
1961
1962 // Distribution before next collision
1963 m_solver->a_oldDistribution(childId, dist) =
1964 (0.5 * (tmpDist1 + tmpDist2) - eqDist_tp1[dist]) * (omega_C / omega_F) * F1B2 + eqDist_tp1[dist];
1965 }
1966
1967 // All cells need to perform collision step next to initialize this variable
1968 m_solver->a_distribution(childId, dist) = m_solver->a_oldDistribution(childId, dist);
1969
1970 sumOfOldDist += m_solver->a_oldDistribution(childId, dist);
1971 }
1972 // Dont forget the rest distribution
1973 m_solver->a_oldDistribution(childId, Ld::lastId()) = rho - sumOfOldDist;
1974 m_solver->a_distribution(childId, Ld::lastId()) = rho - sumOfOldDist;
1975 }
1976 }
1977 }
1978}
1979
1987template <MInt nDim, MInt nDist, class SysEqn>
1988void LbInterfaceDxQy<nDim, nDist, SysEqn>::refineCellCopyPaste(const MInt parentId, const MInt* childIds) {
1989 if(!m_solver->a_isHalo(parentId)) {
1990 for(MInt child = 0; child < IPOW2(nDim); child++) {
1991 if(childIds[child] < 0) continue;
1992 for(MInt v = 0; v < m_solver->noVariables(); v++) {
1993 m_solver->a_variable(childIds[child], v) = m_solver->a_variable(parentId, v);
1994 m_solver->a_oldVariable(childIds[child], v) = m_solver->a_oldVariable(parentId, v);
1995 }
1996 for(MInt dir = 0; dir < nDist; dir++) {
1997 m_solver->a_distribution(childIds[child], dir) = m_solver->a_distribution(parentId, dir);
1998 m_solver->a_oldDistribution(childIds[child], dir) = m_solver->a_oldDistribution(parentId, dir);
1999 }
2000 }
2001 }
2002}
2003
2010template <MInt nDim, MInt nDist, class SysEqn>
2011void LbInterfaceDxQy<nDim, nDist, SysEqn>::removeChildsCopyPaste(const MInt parentId) {
2012 if(!m_solver->a_isHalo(parentId)) {
2013 std::array<MFloat, nDist> forcing{};
2014 getCellForcing(forcing);
2015
2016 // Get children that are to remove
2017 MInt children[IPOW2(nDim)];
2018 for(MInt m = 0; m < IPOW2(nDim); m++) {
2019 children[m] = m_solver->c_childId(parentId, m);
2020 }
2021
2022 m_solver->a_variable(parentId, PV->RHO) = 0.0;
2023 m_solver->a_variable(parentId, PV->U) = 0.0;
2024 m_solver->a_variable(parentId, PV->V) = 0.0;
2025 IF_CONSTEXPR(nDim == 3) m_solver->a_variable(parentId, PV->W) = 0.0;
2026 m_solver->a_oldVariable(parentId, PV->RHO) = 0.0;
2027 m_solver->a_oldVariable(parentId, PV->U) = 0.0;
2028 m_solver->a_oldVariable(parentId, PV->V) = 0.0;
2029 IF_CONSTEXPR(nDim == 3) m_solver->a_oldVariable(parentId, PV->W) = 0.0;
2030 for(MInt i = 0; i < nDist; i++) {
2031 m_solver->a_distribution(parentId, i) = 0.0;
2032 m_solver->a_oldDistribution(parentId, i) = 0.0;
2033 }
2034
2035 // Spatial interpolation: Average all child values
2036 for(MInt m = 0; m < IPOW2(nDim); m++) {
2037 m_solver->a_variable(parentId, PV->RHO) += FFPOW2(nDim) * m_solver->a_variable(children[m], PV->RHO);
2038 m_solver->a_variable(parentId, PV->U) += FFPOW2(nDim) * m_solver->a_variable(children[m], PV->U);
2039 m_solver->a_variable(parentId, PV->V) += FFPOW2(nDim) * m_solver->a_variable(children[m], PV->V);
2040 IF_CONSTEXPR(nDim == 3)
2041 m_solver->a_variable(parentId, PV->W) += FFPOW2(nDim) * m_solver->a_variable(children[m], PV->W);
2042 m_solver->a_oldVariable(parentId, PV->RHO) += FFPOW2(nDim) * m_solver->a_oldVariable(children[m], PV->RHO);
2043 m_solver->a_oldVariable(parentId, PV->U) += FFPOW2(nDim) * m_solver->a_oldVariable(children[m], PV->U);
2044 m_solver->a_oldVariable(parentId, PV->V) += FFPOW2(nDim) * m_solver->a_oldVariable(children[m], PV->V);
2045 IF_CONSTEXPR(nDim == 3)
2046 m_solver->a_oldVariable(parentId, PV->W) += FFPOW2(nDim) * m_solver->a_oldVariable(children[m], PV->W);
2047 for(MInt i = 0; i < nDist; i++) {
2048 m_solver->a_distribution(parentId, i) += FFPOW2(nDim) * m_solver->a_distribution(children[m], i);
2049 m_solver->a_oldDistribution(parentId, i) += FFPOW2(nDim) * m_solver->a_oldDistribution(children[m], i);
2050 }
2051 }
2052 }
2053}
2054
2064template <MInt nDim, MInt nDist, class SysEqn>
2065void LbInterfaceDxQy<nDim, nDist, SysEqn>::removeChildsDupuisFilippova(const MInt parentId) {
2066 if(!m_solver->a_isHalo(parentId)) {
2067 const MInt parentLevel = m_solver->c_level(parentId);
2068 const MInt dParentLevel = m_solver->maxLevel() - parentLevel;
2069
2070 // Relaxation paramter for parent and child refinement level
2071 const MFloat omega_F = F2 / (F1 + 6.0 * m_solver->a_nu(parentId) * FFPOW2(dParentLevel - 1));
2072 const MFloat omega_C = F2 / (F1 + 6.0 * m_solver->a_nu(parentId) * FFPOW2(dParentLevel));
2073 std::array<MFloat, nDist> forcing{};
2074 getCellForcing(forcing);
2075
2076 // Get children that are to remove
2077 MInt children[IPOW2(nDim)];
2078 for(MInt m = 0; m < IPOW2(nDim); m++) {
2079 children[m] = m_solver->c_childId(parentId, m);
2080 }
2081
2082 // Parent and child level have done propagation in the last timestep and will do collision after this...
2083 if((globalTimeStep - 1) % IPOW2(dParentLevel) == 0) {
2084 MFloat rho = F0;
2085 MFloat u[nDim] = {F0};
2086 MFloat rhoOld = F0;
2087 MFloat uOld[nDim] = {F0};
2088 MFloat rho_tp1 = F0;
2089 MFloat u_tp1[nDim] = {F0};
2090 MFloat rho_tm3 = F0;
2091 MFloat u_tm3[nDim] = {F0};
2092
2093 // Spatial interpolation: Sum up all child values
2094 for(MInt m = 0; m < IPOW2(nDim); m++) {
2095 // Macroscopic variables at time t
2096 rho += FFPOW2(nDim) * m_solver->a_variable(children[m], PV->RHO);
2097 for(MInt n = 0; n < nDim; n++) {
2098 u[n] += FFPOW2(nDim) * m_solver->a_variable(children[m], PV->VV[n]);
2099 }
2100 // Macroscopic variables at time t-1
2101 rhoOld += FFPOW2(nDim) * m_solver->a_oldVariable(children[m], PV->RHO);
2102 for(MInt n = 0; n < nDim; n++) {
2103 uOld[n] += FFPOW2(nDim) * m_solver->a_oldVariable(children[m], PV->VV[n]);
2104 }
2105 }
2106
2107 // Temporal extrapolation: f(t+1) = 2*f(t) - f(t-1)
2108 // Needed for initialization of old distribution
2109 rho_tp1 = 2 * rho - rhoOld;
2110 for(MInt d = 0; d < nDim; d++) {
2111 u_tp1[d] = 2 * u[d] - uOld[d];
2112 }
2113
2114 // Temporal extrapolation: f(t-3) = 3*f(t-1) - 2*f(t)
2115 // Needed for initialization of macroscopic variables
2116 rho_tm3 = 3 * rhoOld - 2 * rho;
2117 for(MInt d = 0; d < nDim; d++) {
2118 u_tm3[d] = 3 * uOld[d] - 2 * u[d];
2119 }
2120
2121 // Setting macroscopic variables on parent level
2122 m_solver->a_variable(parentId, PV->RHO) = rhoOld;
2123 for(MInt n = 0; n < nDim; n++) {
2124 m_solver->a_variable(parentId, PV->VV[n]) = uOld[n];
2125 }
2126 m_solver->a_oldVariable(parentId, PV->RHO) = rho_tm3;
2127 for(MInt n = 0; n < nDim; n++) {
2128 m_solver->a_oldVariable(parentId, PV->VV[n]) = u_tm3[n];
2129 }
2130
2131 // Calculate equilibrium for time t+1
2132 std::array<MFloat, nDist> eqDistTp1{};
2133 eqDistTp1 = m_solver->getEqDists(rho_tp1, u_tp1);
2134
2135 // Transform and interpolate post-propagation distributions according to DUPUIS
2136 for(MInt dist = 0; dist < nDist; dist++) {
2137 // Calculate mean distribution from all children
2138 MFloat tmpDist = 0.0;
2139 for(MInt m = 0; m < IPOW2(nDim); m++) {
2140 tmpDist += FFPOW2(nDim) * m_solver->a_oldDistribution(children[m], dist);
2141 }
2142
2143 m_solver->a_oldDistribution(parentId, dist) =
2144 (tmpDist - eqDistTp1[dist]) * (omega_F / omega_C) * 2.0 + eqDistTp1[dist];
2145
2146 // Cell need to perform collision step next to initialize this variable
2147 m_solver->a_distribution(parentId, dist) = m_solver->a_oldDistribution(parentId, dist);
2148 }
2149 } else {
2150 MFloat rho{};
2151 std::array<MFloat, nDim> u{};
2152 MFloat rho_tm1{};
2153 std::array<MFloat, nDim> u_tm1{};
2154
2155 // Spatial interpolation: Sum up all child values
2156 for(MInt m = 0; m < IPOW2(nDim); m++) {
2157 // Macroscopic variables at time t
2158 rho += FFPOW2(nDim) * m_solver->a_variable(children[m], PV->RHO);
2159 for(MInt n = 0; n < nDim; n++) {
2160 u[n] += FFPOW2(nDim) * m_solver->a_variable(children[m], PV->VV[n]);
2161 }
2162 // Macroscopic variables at time t-1
2163 rho_tm1 += FFPOW2(nDim) * m_solver->a_oldVariable(children[m], PV->RHO);
2164 for(MInt n = 0; n < nDim; n++) {
2165 u_tm1[n] += FFPOW2(nDim) * m_solver->a_oldVariable(children[m], PV->VV[n]);
2166 }
2167 }
2168
2169 // Temporal (backward) extrapolation f(t-2) = 2*f(t-1) - f(t)
2170 MFloat rho_tm2{};
2171 std::array<MFloat, nDim> u_tm2{};
2172 rho_tm2 = 2 * rho_tm1 - 1 * rho;
2173 for(MInt d = 0; d < nDim; d++) {
2174 u_tm2[d] = 2 * u_tm1[d] - 1 * u[d];
2175 }
2176
2177 // Setting macroscopic variables for parent cell
2178 m_solver->a_variable(parentId, PV->RHO) = rho;
2179 for(MInt d = 0; d < nDim; d++) {
2180 m_solver->a_variable(parentId, PV->VV[d]) = u[d];
2181 }
2182 m_solver->a_oldVariable(parentId, PV->RHO) = rho_tm2;
2183 for(MInt d = 0; d < nDim; d++) {
2184 m_solver->a_oldVariable(parentId, PV->VV[d]) = u_tm2[d];
2185 }
2186
2187 // Calculate equilibirum for time t
2188 std::array<MFloat, nDist> eqDist{};
2189 eqDist = m_solver->getEqDists(rho, u.data());
2190
2191 // Transform and interpolate post-collision distributions according to FILIPPOVA
2192 const MFloat tau_F = 1 / omega_F;
2193 const MFloat tau_C = 1 / omega_C;
2194 for(MInt dist = 0; dist < nDist; dist++) {
2195 // Calculate mean distribution from all children
2196 MFloat tmpDist = 0.0;
2197 for(MInt m = 0; m < IPOW2(nDim); m++) {
2198 tmpDist += FFPOW2(nDim) * m_solver->a_distribution(children[m], dist);
2199 }
2200
2201 m_solver->a_distribution(parentId, dist) =
2202 (tmpDist - eqDist[dist]) * ((tau_C - 1) / (tau_F - 1)) * 2.0 + eqDist[dist];
2203
2204 m_solver->a_oldDistribution(parentId, dist) = m_solver->a_distribution(parentId, dist);
2205 }
2206 }
2207 }
2208}
LbInterfaceDxQy::restriction10
virtual void restriction10()
Definition: lbinterfacedxqy.cpp:602
LbInterfaceDxQy::removeChildsDupuisFilippova
virtual void removeChildsDupuisFilippova(const MInt parentId)
Initialize parent variables from children.
Definition: lbinterfacedxqy.cpp:2065
LbInterfaceDxQy::setInterfaceFunctions
void setInterfaceFunctions()
Setting function pointer to the chosen interface treatment method.
Definition: lbinterfacedxqy.cpp:78
LbInterfaceDxQy::restrictionThermalRohde
virtual void restrictionThermalRohde()
Fine to coarse grid fot thermal LB.
Definition: lbinterfacedxqy.cpp:1622
LbInterfaceDxQy::LbInterfaceDxQy
LbInterfaceDxQy(LbSolver< nDim > *solver)
LbInterfaceDxQy::prolongationThermalRohde
virtual void prolongationThermalRohde()
Coarse to fine grid for thermal LB.
Definition: lbinterfacedxqy.cpp:1494
LbInterfaceDxQy::restrictionDupuis
virtual void restrictionDupuis()
Definition: lbinterfacedxqy.cpp:1310
LbInterfaceDxQy::prolongationRohde
virtual void prolongationRohde()
Definition: lbinterfacedxqy.cpp:1416
LbInterfaceDxQy::prolongationDupuisCompressible
void prolongationDupuisCompressible()
Definition: lbinterfacedxqy.cpp:965
LbInterfaceDxQy::prolongation10
virtual void prolongation10()
Definition: lbinterfacedxqy.cpp:143
LbInterfaceDxQy::prolongationDupuis
virtual void prolongationDupuis()
Definition: lbinterfacedxqy.cpp:960
LbInterfaceDxQy::setAdaptationFunctions
void setAdaptationFunctions()
Setting the adaptation functions chosen via property.
Definition: lbinterfacedxqy.cpp:1670
LbInterfaceDxQy::refineCellDupuis
virtual void refineCellDupuis(const MInt parentId, const MInt *childIds)
Initialize child variables from parent.
Definition: lbinterfacedxqy.cpp:1714
LbInterfaceDxQy::restrictionRohde
virtual void restrictionRohde()
Definition: lbinterfacedxqy.cpp:1576
LbInterfaceDxQy::prolongationThermalDupuis
virtual void prolongationThermalDupuis()
Coarse to fine grid for thermal LB.
Definition: lbinterfacedxqy.cpp:978
LbInterfaceDxQy::removeChildren
virtual void removeChildren(const MInt parentId) override final
Definition: lbinterfacedxqy.cpp:1700
LbInterfaceDxQy::~LbInterfaceDxQy
virtual ~LbInterfaceDxQy()
D'tor for the interface class.
LbInterfaceDxQy::restrictionThermalDupuis
virtual void restrictionThermalDupuis()
Fine to coarse grid for thermal LB.
Definition: lbinterfacedxqy.cpp:1328
LbInterfaceDxQy::refineCellCopyPaste
virtual void refineCellCopyPaste(const MInt parentId, const MInt *childIds)
Initialize child variables from parent.
Definition: lbinterfacedxqy.cpp:1988
LbInterfaceDxQy::prolongationDupuis_
void prolongationDupuis_()
Coarse to fine.
Definition: lbinterfacedxqy.cpp:772
LbInterfaceDxQy::prolongation
void prolongation()
Performing the chosen prolongation method.
Definition: lbinterfacedxqy.cpp:122
LbInterfaceDxQy::refineCell
virtual void refineCell(const MInt parentId, const MInt *childIds) override final
Definition: lbinterfacedxqy.cpp:1695
LbInterfaceDxQy::restrictionDupuisCompressible
void restrictionDupuisCompressible()
Definition: lbinterfacedxqy.cpp:1315
LbInterfaceDxQy::restrictionDupuis_
void restrictionDupuis_()
Fine to coarse grid.
Definition: lbinterfacedxqy.cpp:1235
LbInterfaceDxQy::restriction
void restriction()
Performing the chosen restriction method.
Definition: lbinterfacedxqy.cpp:132
LbInterfaceDxQy::removeChildsCopyPaste
virtual void removeChildsCopyPaste(const MInt parentId)
Initialize parent variables from children.
Definition: lbinterfacedxqy.cpp:2011
LbSolverDxQy
This class represents all LB models.
Definition: lbsolverdxqy.h:29
string2enum
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
ROHDE
@ ROHDE
Definition: enums.h:258
FILIPPOVA
@ FILIPPOVA
Definition: enums.h:258
INIT_COPYPASTE
@ INIT_COPYPASTE
Definition: enums.h:261
INIT_DUPUIS_FILIPPOVA
@ INIT_DUPUIS_FILIPPOVA
Definition: enums.h:261
mTerm
void mTerm(const MInt errorCode, const MString &location, const MString &message)
Definition: functions.cpp:29
globalTimeStep
MInt globalTimeStep
Definition: globalvariables.cpp:33
LbCell::WasNewlyCreated
@ WasNewlyCreated
IPOW2
constexpr MLong IPOW2(MInt x)
Definition: maiaconstants.h:154
FPOW2
constexpr MFloat FPOW2(MInt x)
Definition: maiaconstants.h:158
FFPOW2
constexpr MFloat FFPOW2(MInt x)
Definition: maiaconstants.h:162
MInt
int32_t MInt
Definition: maiatypes.h:62
MFloat
double MFloat
Definition: maiatypes.h:52
MBool
bool MBool
Definition: maiatypes.h:58
id
MInt id
Definition: maiatypes.h:71
dist
MFloat dist(const Point< DIM > &p, const Point< DIM > &q)
Definition: pointbox.h:54
a
Definition: contexttypes.h:19