MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
couplingdgape.h
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#ifndef COUPLINGDGAPE_H_
8#define COUPLINGDGAPE_H_
9
10#include "DG/dgcartesiangalerkinprojection.h"
11#include "DG/dgcartesiansolver.h"
12#include "UTIL/parallelfor.h"
13#include "coupling.h"
14
15/* idea
16 +----------+
17 +---------->| Coupling |<---------+
18 | +----+-----+ |
19 | ^ |
20 +------+-----+ +-----+------+ +-----+------+
21 | CouplingFv | | CouplingDg | | CouplingLb |
22 +------+-----+ +-----+------+ +-----+------+
23 ^ ^ ^
24 | | |
25 | +------+--------+ |
26 | +--->| CouplingDgApe |<--+ |
27 | | +---------------+ | |
28 | | | |
29 +--------+----+-------+ +-+----+----------+
30 | DgCcAcousticPerturb | | CouplingLbDgApe |
31 +---------------------+ +-----------------+
32*/
33
34template <MInt nDim_>
35struct MeanVariables;
36template <>
37struct MeanVariables<2> {
38 static constexpr const MInt nDim = 2;
39 static constexpr const MInt noVorticities = 1;
40
41 // Mean Lamb vector
42 static constexpr const MInt LAMB0_X = 0;
43 static constexpr const MInt LAMB0_Y = 1;
44 static constexpr const MInt LAMB0[nDim] = {0, 1};
45
46 // Mean velocities
47 static constexpr const MInt U0 = 2;
48 static constexpr const MInt V0 = 3;
49 static constexpr const MInt UU0[nDim] = {2, 3};
50
51 static constexpr const MInt RHO0 = 4;
52 static constexpr const MInt P0 = 5;
53 static constexpr const MInt C0 = 6;
54
55 // Mean speed of sound
56 static constexpr const MInt DC0_X = 7;
57 static constexpr const MInt DC0_Y = 8;
58 static constexpr const MInt DC0[nDim] = {7, 8};
59
60 // Mean vorticity
61 static constexpr const MInt VORT0_Z = 9;
62 static constexpr const MInt VORT0[noVorticities] = {9};
63
64 // Mean velocity gradient
65 static constexpr const MInt DU_DX = 10;
66 static constexpr const MInt DU_DY = 11;
67
68 static constexpr const MInt DV_DX = 12;
69 static constexpr const MInt DV_DY = 13;
70
71 static constexpr const MInt GRADU[nDim * nDim] = {10, 11, 12, 13};
72
73 // du/dx, dv/dy, dw/dz for the divergence
74 // NOTE: be careful when accessing data since entries of DU are not consecutive in memory if the
75 // full mean velocity gradient is loaded
76 static constexpr const MInt DU[nDim] = {10, 13};
77
78 // Mean gradient of rho
79 static constexpr const MInt DRHO_DX = 14;
80 static constexpr const MInt DRHO_DY = 15;
81 static constexpr const MInt DRHO[nDim] = {14, 15};
82
83 // Mean gradient of p
84 static constexpr const MInt DP_DX = 16;
85 static constexpr const MInt DP_DY = 17;
86 static constexpr const MInt DP[nDim] = {16, 17};
87
88 // Mean gradient of rho*div(u)
89 static constexpr const MInt RHODIVU_X = 18;
90 static constexpr const MInt RHODIVU_Y = 19;
91 static constexpr const MInt RHODIVU[nDim] = {18, 19};
92
93 // Mean gradient of u*grad(rho)
94 static constexpr const MInt UGRADRHO_X = 20;
95 static constexpr const MInt UGRADRHO_Y = 21;
96 static constexpr const MInt UGRADRHO[nDim] = {20, 21};
97
98 // Mean of (gradient of p divided by rho)
99 static constexpr const MInt GRADPRHO_X = 22;
100 static constexpr const MInt GRADPRHO_Y = 23;
101 static constexpr const MInt GRADPRHO[nDim] = {22, 23};
102
103 // Sum of products of velocity components with corresponding velocity gradients
104 static constexpr const MInt UGRADU_X = 24;
105 static constexpr const MInt UGRADU_Y = 25;
106 static constexpr const MInt UGRADU[nDim] = {24, 25};
107
108 // Note: If the totalNoMeanVars value is changed here, it must also be changed
109 // for s_totalNoMeanVars
110 static constexpr const MInt totalNoMeanVars = 26;
111 static const std::array<MString, totalNoMeanVars> names;
112 static const std::array<MString, 6> nodeVarNames;
113};
114template <>
115struct MeanVariables<3> {
116 static constexpr const MInt nDim = 3;
117 static constexpr const MInt noVorticities = 3;
118
119 // Mean lamb vector
120 static constexpr const MInt LAMB0_X = 0;
121 static constexpr const MInt LAMB0_Y = 1;
122 static constexpr const MInt LAMB0_Z = 2;
123 static constexpr const MInt LAMB0[nDim] = {0, 1, 2};
124 // Mean velocities
125 static constexpr const MInt U0 = 3;
126 static constexpr const MInt V0 = 4;
127 static constexpr const MInt W0 = 5;
128 static constexpr const MInt UU0[nDim] = {3, 4, 5};
129
130 static constexpr const MInt RHO0 = 6;
131 static constexpr const MInt P0 = 7;
132 static constexpr const MInt C0 = 8;
133
134 // Mean speed of sound
135 static constexpr const MInt DC0_X = 9;
136 static constexpr const MInt DC0_Y = 10;
137 static constexpr const MInt DC0_Z = 11;
138 static constexpr const MInt DC0[nDim] = {9, 10, 11};
139
140 // Mean vorticities
141 static constexpr const MInt VORT0_X = 12;
142 static constexpr const MInt VORT0_Y = 13;
143 static constexpr const MInt VORT0_Z = 14;
144 static constexpr const MInt VORT0[noVorticities] = {12, 13, 14};
145
146 // Mean velocity gradient
147 static constexpr const MInt DU_DX = 15;
148 static constexpr const MInt DU_DY = 16;
149 static constexpr const MInt DU_DZ = 17;
150
151 static constexpr const MInt DV_DX = 18;
152 static constexpr const MInt DV_DY = 19;
153 static constexpr const MInt DV_DZ = 20;
154
155 static constexpr const MInt DW_DX = 21;
156 static constexpr const MInt DW_DY = 22;
157 static constexpr const MInt DW_DZ = 23;
158
159 static constexpr const MInt GRADU[nDim * nDim] = {15, 16, 17, 18, 19, 20, 21, 22, 23};
160
161 // du/dx, du/dy, dw/dz for the divergence
162 // NOTE: be careful when accessing data since entries of DU are not consecutive in memory if the
163 // full mean velocity gradient is loaded
164 static constexpr const MInt DU[nDim] = {15, 19, 23};
165
166 // Mean gradient of rho
167 static constexpr const MInt DRHO_DX = 24;
168 static constexpr const MInt DRHO_DY = 25;
169 static constexpr const MInt DRHO_DZ = 26;
170 static constexpr const MInt DRHO[nDim] = {24, 25, 26};
171
172 // Mean gradient of p
173 static constexpr const MInt DP_DX = 27;
174 static constexpr const MInt DP_DY = 28;
175 static constexpr const MInt DP_DZ = 29;
176 static constexpr const MInt DP[nDim] = {27, 28, 29};
177
178 // Mean gradient of rho*div(u)
179 static constexpr const MInt RHODIVU_X = 30;
180 static constexpr const MInt RHODIVU_Y = 31;
181 static constexpr const MInt RHODIVU_Z = 32;
182 static constexpr const MInt RHODIVU[nDim] = {30, 31, 32};
183
184 // Mean gradient of u*grad(rho)
185 static constexpr const MInt UGRADRHO_X = 33;
186 static constexpr const MInt UGRADRHO_Y = 34;
187 static constexpr const MInt UGRADRHO_Z = 35;
188 static constexpr const MInt UGRADRHO[nDim] = {33, 34, 35};
189
190 // Mean of (gradient of p divided by rho)
191 static constexpr const MInt GRADPRHO_X = 36;
192 static constexpr const MInt GRADPRHO_Y = 37;
193 static constexpr const MInt GRADPRHO_Z = 38;
194 static constexpr const MInt GRADPRHO[nDim] = {36, 37, 38};
195
196 // Sum of products of velocity components with corresponding velocity gradients
197 static constexpr const MInt UGRADU_X = 39;
198 static constexpr const MInt UGRADU_Y = 40;
199 static constexpr const MInt UGRADU_Z = 41;
200 static constexpr const MInt UGRADU[nDim] = {39, 40, 41};
201
202 // Note: If the totalNoMeanVars value is changed here, it must also be changed
203 // for s_totalNoMeanVars
204 static constexpr const MInt totalNoMeanVars = 42;
205 static const std::array<MString, totalNoMeanVars> names;
206 static const std::array<MString, 8> nodeVarNames;
207};
208
213template <MInt nDim, class CouplingDonor>
214class CouplingDgApe : public CouplingDonor, public CouplingDg<nDim, DgSysEqnAcousticPerturb<nDim>> {
215 // Typedefs
216 public:
217 using Base = Coupling;
218 using Base::startLoadTimer;
219 using Base::stopLoadTimer;
220
221 using BaseDonor = CouplingDonor;
222 using DonorSolverType = typename BaseDonor::solverType;
223
224 using SysEqn = DgSysEqnAcousticPerturb<nDim>;
225 using ProjectionType = DgGalerkinProjection<nDim>;
226 using BaseDg = CouplingDg<nDim, DgSysEqnAcousticPerturb<nDim>>;
227 using DgCartesianSolverType = typename BaseDg::solverType;
228
229 using BaseDg::dgSolver;
230 using BaseDg::elements;
231 using BaseDg::externalSource;
232 using BaseDg::maxPolyDeg;
233 using BaseDg::minPolyDeg;
234 using BaseDg::noElements;
235 using BaseDg::returnIdleRecord;
236 using BaseDg::returnLoadRecord;
237 using BaseDg::saveNodalData;
238 using BaseDg::solverId; // TODO labels:COUPLER,DG use couplingId instead of solverId to request coupling related
239 // properties?
240 using BaseDg::sysEqn;
241 using MV = MeanVariables<nDim>;
242
243 //--Ctor & Dtor
244 CouplingDgApe(const MInt couplingId, DgCartesianSolverType* dg, DonorSolverType* ds);
245 virtual ~CouplingDgApe() = default;
246
247 //--Methods
248 // Virtual methods to override
249 void init() override {
250 readProperties();
251 initCoupler();
252 initData();
253 };
254 void finalizeSubCoupleInit(MInt /*couplingStep*/) override final{};
255 // void finalizeCouplerInit() override{}; // defined in child classes
256 void preCouple(MInt /*recipestep*/) final;
257 void subCouple(MInt /*recipeStep*/, MInt /*solverId*/, std::vector<MBool>& /*solverCompleted*/) override final{};
258 void postCouple(MInt /*recipestep*/) override final{};
259 void cleanUp() override final{};
260
261 void getCouplingTimings(std::vector<std::pair<MString, MFloat>>& timings, const MBool allTimings) override {
262 const MString namePrefix = "c" + std::to_string(this->couplerId()) + "_";
263
264 timings.emplace_back(namePrefix + "loadCouplerApe", returnLoadRecord());
265 timings.emplace_back(namePrefix + "idleCouplerApe", returnIdleRecord());
266
267#ifdef MAIA_TIMER_FUNCTION
268 timings.emplace_back(namePrefix + "preCouple", RETURN_TIMER_TIME(m_timers[Timers::PreCouple]));
269
270 if(allTimings) {
271 // Full set of timings
272 timings.emplace_back(namePrefix + "storeSources", RETURN_TIMER_TIME(m_timers[Timers::StoreSources]));
273 timings.emplace_back(namePrefix + "calcSources", RETURN_TIMER_TIME(m_timers[Timers::CalcSources]));
274 timings.emplace_back(namePrefix + "projectSourceTerms", RETURN_TIMER_TIME(m_timers[Timers::ProjectSourceTerms]));
275
276 timings.emplace_back(namePrefix + "applyStoredSources", RETURN_TIMER_TIME(m_timers[Timers::ApplyStoredSources]));
277 } else {
278 // Reduced/essential set of timings
279 timings.emplace_back(namePrefix + "projectSourceTerms", RETURN_TIMER_TIME(m_timers[Timers::ProjectSourceTerms]));
280 }
281#endif
282 };
283
284 void getDomainDecompositionInformation(std::vector<std::pair<MString, MInt>>& domainInfo) override final {
285 const MString namePrefix = "c" + std::to_string(this->couplerId()) + "_";
286
287 const MInt noCoupledDgElements = m_calcSourceElements.size();
288 const MInt noCoupledDonorCells = m_calcSourceDonorCells.size();
289
290 domainInfo.emplace_back(namePrefix + "noCoupledDgElements", noCoupledDgElements);
291 domainInfo.emplace_back(namePrefix + "noCoupledDonorCells", noCoupledDonorCells);
292 };
293 MInt noCouplingTimers(const MBool allTimings) const override {
294#ifdef MAIA_TIMER_FUNCTION
295 if(allTimings) {
296 return 7;
297 } else {
298 return 4;
299 }
300#else
301 return 2;
302#endif
303 };
304
305 //
306 constexpr static MInt noVelocities() { return nDim; }
307 constexpr static MInt noVorticities() { return (nDim == 2) ? 1 : 3; }
308 MInt noMeanVars() const { return m_activeMeanVars.size(); }
309
310 protected:
311 //--Methods
312 void calcInitialCondition(const MFloat time);
313 void initCoupler();
314 void initData();
315 void initProjection();
316 void initRestart(const MFloat time, const MFloat dt);
317 void initSourceFilter();
318 void initTimers();
319
320 MBool loadCouplingData(const MString& filename, const MString& name_, const MInt stride, MFloat* data);
321 void loadMeanQuantities(const MBool skipNodeVars = false);
322 void projectToElement(const MInt elementId, const MInt noVars, const MFloat* data, const MFloat* defaultValues,
323 MFloat* target);
324 void saveFilterValues();
325 void saveFilterValuesDonor();
326
327 // Source term calculation
328 void applyStoredSources(const MFloat time);
329 void storeSources(const MFloat time, const MInt timeStep);
330
331 virtual void performUnitConversion(const MString& /*name*/, const MInt /*count*/, const MInt /*stride*/,
332 MFloat* /*data*/){};
333 virtual void calcSourceLambLinearized(const MFloat* const velocity, const MFloat* const vorticity,
334 MFloat* sourceTerms) = 0;
335 virtual void calcSourceLambNonlinear(const MFloat* const velocity, const MFloat* const vorticity,
336 MFloat* const sourceTerms) = 0;
337 virtual void calcSourceQmII(const MFloat* const velocity, MFloat* const sourceTerms) = 0;
338 virtual void calcSourceQmIII(const MFloat* const velocity, MFloat* sourceTerms) = 0;
339 virtual void calcSourceQe(const MFloat* const velocity, const MFloat time, MFloat* const sourceTerms) = 0;
340 virtual void calcSourceQc(const MFloat* const velocity, MFloat* const sourceTerms, const MFloat time,
341 const MInt timeStep) = 0;
342
343 void saveSourceTermsDonorGrid(const MInt timeStep, const MFloat* const data);
344 void saveSourceTermsTargetGrid(const MInt timeStep);
345
346 // Virtual methods to be overriden
347 void readProperties() override final;
348 void checkProperties() override final{};
349
350 // Donor solver specific 'accessors'
351 virtual DonorSolverType& donorSolver(const MInt solverId = 0) const = 0;
352 virtual void getDonorVelocityAndVorticity(const std::vector<MInt>& donorCellIds, MFloatScratchSpace& p_velocity,
353 MFloatScratchSpace& p_vorticity) = 0;
354
355 //--Variables
356 MBool m_isRestart = false;
357 MBool m_hasDgCartesianSolver = false;
358 MBool m_hasDonorCartesianSolver = false;
359
360 MInt m_maxNoNodesXD = -1;
361
362 // Time step calculation
363 // MInt m_timeStepOffset = -1; ///< Time step offset between CAA simulation and LES simulation
364 MFloat m_fixedTimeStep = -1.0;
365 MFloat m_previousTime =
366 std::numeric_limits<MFloat>::quiet_NaN();
367 MInt m_previousTimeStep = -1;
368
369 std::vector<MInt> m_activeMeanVars{};
370 MString m_meanDataFileName{};
371 std::vector<MInt> m_activeSourceTerms{};
372 MBool m_saveSourceTermsDonorGrid = false;
373 MInt m_saveSourceTermsInterval = -1;
374 MBool m_calcProjectionError = false;
375
376 // Donor projection information
377 MInt m_noActiveDonorCells = -1;
379 MInt m_noDonorCells = -1;
380 std::vector<std::vector<MInt>> m_elementDonorLevels = {};
381 std::vector<std::vector<MInt>> m_elementDonorMap = {};
382 std::vector<std::vector<MInt>> m_elementDonorPos = {};
383
384 // Projection
385 std::vector<ProjectionType> m_projection{};
386
387 // Source term calculation
388 MBool m_isFirstSourceCalculation = true;
389 std::vector<MFloat> m_sourceTerms{};
390 std::vector<MInt> m_calcSourceDonorCells{};
391 std::vector<MInt> m_calcSourceElements{};
392 MFloat m_sourceFactor = 1.0;
393
394 // Error checking
395 MBool m_checkConservation = false;
396 MFloat m_maxConservationError = -1.0;
397 MFloat m_maxL2Error = -1.0;
398
399 // Mean variables
400 std::vector<MFloat> m_meanVars{};
401
402 // Geometry
403 std::vector<MInt> m_elementInsideBody{};
404
405 // Filter
406 // MBool m_applySourceTermFilter = true; ///< Store whether to apply filtering to calculated source terms
407 Filter<nDim> m_sourceTermFilter;
408 std::vector<MFloat> m_filter{};
409 std::vector<MFloat> m_filterDonor{};
410 MInt m_noCutModesLowPass = 0;
411 MBool m_useLowPassFilter = false;
412 std::vector<MFloatTensor> m_vdmLowPass{};
413 MBool m_applySourceFilterDonor = false;
414 MBool m_saveSourceTermFilter = false;
415
416 MBool m_projectionFilter = false;
417 std::vector<MFloat> m_projectionFilterBox{};
418
419 //--Structs
420 // Total number of mean variables (Note: MV::totalNoMeanVars can't be used).
421 // For the list of mean variables, see class definition for MV below.
422 static const MInt s_totalNoMeanVars = MV::totalNoMeanVars;
423 // List of indices of all mean variables in m_activeMeanVars, i.e. for each
424 // active mean variable X, its position in m_activeMeanVars is given by
425 // m_meanVarsIndex[MV::X]
426 std::array<MInt, s_totalNoMeanVars> m_meanVarsIndex{};
427 // static const std::array<MString, s_totalNoMeanVars> s_meanVarNames; ///< Names for all possible mean variables in
428 // MV
429 struct ST;
430 static const std::array<MString, ST::totalNoSourceTerms> s_sourceTermNames;
431
432 // Timing
433 struct Timers {
434 enum {
435 // Timer group and main timer
436 TimerGroup,
437 Class,
438
439 // Regular (method-specific timers)
440 Constructor,
441 Init,
442 InitialCondition,
443 LoadMeanQuantities,
444 PreCouple,
445 StoreSources,
446 LoadInstantaneous,
447 CalcSources,
448 CalcSourceLamb,
449 CalcSourceQmII,
450 CalcSourceQmIII,
451 CalcSourceQe,
452 CalcSourceQc,
453 ProjectSourceTerms,
454 ApplyStoredSources,
455
456 // Accumulated timers for certain subsystems
457 Accumulated,
458 LoadCouplingData,
459
460 // Counter
461 _count
462 };
463 };
464 std::array<MInt, Timers::_count> m_timers;
465
466 private:
467 void neededMeanVarsForSourceTerm(const MInt sourceTerm, std::vector<MInt>& meanVars) const;
468};
469
470// Use to differentiate between different source terms
471template <MInt nDim, class CouplingDonor>
472struct CouplingDgApe<nDim, CouplingDonor>::ST {
473 static constexpr const MInt Q_mI = 0;
474 static constexpr const MInt Q_mI_linear = 1;
475 static constexpr const MInt Q_mII = 2;
476 static constexpr const MInt Q_mIII = 3;
477 static constexpr const MInt Q_mIV = 4;
478 static constexpr const MInt Q_c = 5;
479 static constexpr const MInt Q_e = 6;
480
481 static constexpr const MInt totalNoSourceTerms = 7;
482};
483
484// Source term names
485template <MInt nDim, class CouplingDonor>
486const std::array<MString, CouplingDgApe<nDim, CouplingDonor>::ST::totalNoSourceTerms>
487 CouplingDgApe<nDim, CouplingDonor>::s_sourceTermNames = {
488 {"q_mI", "q_mI_linear", "q_mII", "q_mIII", "q_mIV", "q_c", "q_e"}};
489
490//--Ctor & Dtor-----------------------------------------------------------------
491
499template <MInt nDim, class CouplingDonor>
500CouplingDgApe<nDim, CouplingDonor>::CouplingDgApe(const MInt couplingId, DgCartesianSolverType* dg, DonorSolverType* ds)
501 : Coupling(couplingId),
502 BaseDonor(couplingId, ds),
503 CouplingDg<nDim, SysEqn>(couplingId, dg),
504 m_sourceTermFilter(dg->solverId()) {
505 TRACE();
506
507 // Ensure consistency in arguments
508 static_assert(MV::totalNoMeanVars == s_totalNoMeanVars, "Total number of mean vars is inconsistent");
509
510 initTimers();
511
512 RECORD_TIMER_STOP(m_timers[Timers::Constructor]);
513}
514
515//--Methods (public)------------------------------------------------------------
516
518template <MInt nDim, class CouplingDonor>
519void CouplingDgApe<nDim, CouplingDonor>::preCouple(MInt) {
520 TRACE();
521 RECORD_TIMER_START(m_timers[Timers::PreCouple]);
522
523 const MInt timeStep = (m_hasDgCartesianSolver) ? dgSolver().getCurrentTimeStep() : -1;
524 const MFloat time = (m_hasDgCartesianSolver) ? dgSolver().time() : -1.0;
525
526 if(m_hasDgCartesianSolver && timeStep == m_previousTimeStep) {
527 TERMM(1, "Error: preCouple already called for this time step.");
528 }
529
530 // Calculate source terms
531 storeSources(time, timeStep);
532 // Add stored source terms to the external source buffer of the DG solver
533 applyStoredSources(time);
534 // Update time step
535 m_previousTimeStep = timeStep;
536
537 RECORD_TIMER_STOP(m_timers[Timers::PreCouple]);
538}
539
540//--Methods (protected)---------------------------------------------------------
541
548template <MInt nDim, class CouplingDonor>
549void CouplingDgApe<nDim, CouplingDonor>::calcInitialCondition(const MFloat NotUsed(time)) {
550 TRACE();
551 RECORD_TIMER_START(m_timers[Timers::InitialCondition]);
552
553 // Load mean quantities from disk into node variables
554 loadMeanQuantities();
555
556 RECORD_TIMER_STOP(m_timers[Timers::InitialCondition]);
557}
558
565template <MInt nDim, class CouplingDonor>
566void CouplingDgApe<nDim, CouplingDonor>::initCoupler() {
567 TRACE();
568 RECORD_TIMER_START(m_timers[Timers::Init]);
569
570 m_hasDgCartesianSolver = dgSolver().isActive();
571 if(!m_hasDgCartesianSolver) {
572 TERMM_IF_NOT_COND(noElements() == 0, "inactive DG solver should have no elements");
573 }
574
575 // Determine the total number of donor cells on this domain, i.e., total number of donor cells on
576 // this domain
577 m_noDonorCells = donorSolver().noInternalCells();
578 if(m_noDonorCells < 0) { // Reset if donorSolver is inactive
579 m_noDonorCells = 0;
580 }
581 m_hasDonorCartesianSolver = (m_noDonorCells > 0);
582 TERMM_IF_NOT_COND(m_hasDonorCartesianSolver == donorSolver().isActive(), "Donor solver status error");
583
584 // Calculate maximum number of nodes
585 m_maxNoNodesXD = (m_hasDgCartesianSolver) ? ipow(maxPolyDeg() + 1, nDim) : -1;
586
587 // Initialize the Galerkin projection(s) and determine the information needed
588 // to perform the Galerkin projection for all elements that have donor cells
589 initProjection();
590
591 // Allocate persistent storage for the source filter values. The maximum
592 // polynomial degree is used here to be consisent with the memory allocation
593 // in the solver.
594 m_filter.resize(std::max(noElements() * m_maxNoNodesXD, 0));
595 // Storage for donor cell filter
596 m_filterDonor.resize(m_noDonorCells);
597
598 // Allocate storage for element ids for which coupling source terms need to be
599 // computed
600 m_calcSourceElements.reserve(noElements());
601
602 // Calculate source filter values and store element ids on which source terms
603 // need to be computed, i.e. elements that have nonzero filter value(s) AND
604 // also have donor cell(s)
605 initSourceFilter();
606
607 // Save filter values to file
608 if(m_saveSourceTermFilter) {
609 saveFilterValues();
610 saveFilterValuesDonor();
611 }
612
613 // Allocate persistent storage for the source terms of all elements for which
614 // coupling source terms need to be computed. The maximum polynomial degree is
615 // used here to be consistent with the memory allocation in the solver.
616 m_sourceTerms.resize(std::max(static_cast<MInt>(m_calcSourceElements.size()) * m_maxNoNodesXD * SysEqn::noVars(), 0));
617
618 // Determine the set of donor cells on which source terms need to be computed
619 std::set<MInt> calcSourceDonorCells;
620 for(MInt elementId : m_calcSourceElements) {
621 calcSourceDonorCells.insert(m_elementDonorMap[elementId].begin(), m_elementDonorMap[elementId].end());
622 }
623 // Number of 'active' donor cells
624 m_noActiveDonorCells = calcSourceDonorCells.size();
625 if(!m_hasDgCartesianSolver) {
626 TERMM_IF_COND(m_noActiveDonorCells > 0, "Inactive DG solver cannot have active donor cells.");
627 }
628
629 // Store 'active' donor cell ids in m_calcSourceDonorCells
630 m_calcSourceDonorCells.resize(m_noActiveDonorCells);
631 m_calcSourceDonorCells.assign(calcSourceDonorCells.begin(), calcSourceDonorCells.end());
632
633 // Allocate persistent storage for the mean variables of the donor cells
634 m_meanVars.resize(m_noActiveDonorCells * noMeanVars());
635
636 if(m_noCutModesLowPass > 0 && m_hasDgCartesianSolver) {
638 // Temporarily create interpolation objects
640
641 // Check for SBP Mode
642 const MInt defaultSbpMode = 0;
643 const MInt sbpMode = Context::getSolverProperty("sbpMode", solverId(), AT_, &defaultSbpMode);
644
645 // Determine SBP Operator
646 const MString defaultSbpOperator = "";
647 const MString sbpOperator =
648 Context::getSolverProperty<MString>("sbpOperator", solverId(), AT_, &defaultSbpOperator);
649
650 // Determine initial number of nodes
651 MInt initNoNodes1D = -1;
652 initNoNodes1D = Context::getSolverProperty<MInt>("initNoNodes", solverId(), AT_, &initNoNodes1D);
653
654 // Determine DG integration method (same as in dgsolver.cpp)
655 const MString defaultDgIntegrationMethod = "DG_INTEGRATE_GAUSS";
656 const MInt dgIntegrationMethod = string2enum(
657 Context::getSolverProperty<MString>("dgIntegrationMethod", solverId(), AT_, &defaultDgIntegrationMethod));
658
659 // Determine DG polynomial type (same as in dgsolver.cpp)
660 const MString defaultDgPolynomialType = "DG_POLY_LEGENDRE";
661 const MInt dgPolynomialType =
662 string2enum(Context::getSolverProperty<MString>("dgPolynomialType", solverId(), AT_, &defaultDgPolynomialType));
663
664 // Convert integers to enums
665 auto polyType = static_cast<DgPolynomialType>(dgPolynomialType);
666 auto intMethod = static_cast<DgIntegrationMethod>(dgIntegrationMethod);
667
668 // Create & init interpolation objects
669 std::vector<DgInterpolation> interpolation(maxPolyDeg() + 1);
670 for(MInt i = std::max(0, minPolyDeg() - m_noCutModesLowPass); i < maxPolyDeg() + 1; i++) {
671 const MInt noNodes1D = sbpMode ? initNoNodes1D : (i + 1);
672 interpolation[i].init(i, polyType, noNodes1D, intMethod, sbpMode, sbpOperator);
673 }
674
676 // Create & store Vandermonde matrices
678
679 // Create
680 m_vdmLowPass.clear();
681 m_vdmLowPass.resize(maxPolyDeg() + 1 - m_noCutModesLowPass);
682
683 // Create matrices
684 for(MInt polyDeg = std::max(0, minPolyDeg() - m_noCutModesLowPass);
685 polyDeg < maxPolyDeg() + 1 - m_noCutModesLowPass;
686 polyDeg++) {
687 const MInt noNodesIn = polyDeg + 1;
688 const MInt noNodesOut = polyDeg + 1 + m_noCutModesLowPass;
689 m_vdmLowPass[polyDeg].resize(noNodesOut, noNodesIn);
690 const DgInterpolation& interpIn = interpolation[polyDeg];
691 const DgInterpolation& interpOut = interpolation[polyDeg + m_noCutModesLowPass];
692 maia::dg::interpolation::calcPolynomialInterpolationMatrix(noNodesIn,
693 &interpIn.m_nodes[0],
694 noNodesOut,
695 &interpOut.m_nodes[0],
696 &interpIn.m_wBary[0],
697 &m_vdmLowPass[polyDeg][0]);
698 }
699 } else {
700 m_vdmLowPass.clear();
701 }
702
703 RECORD_TIMER_STOP(m_timers[Timers::Init]);
704}
705
709template <MInt nDim, class CouplingDonor>
710void CouplingDgApe<nDim, CouplingDonor>::initData() {
711 TRACE();
712 RECORD_TIMER_START(m_timers[Timers::Init]);
713
714 if(dgSolver().m_restart) {
715 // Load mean vars for source term calculation at a restart
716 const MFloat time = dgSolver().time();
717 // TODO labels:COUPLER,DG,toremove @ansgar_unified is dt already set?
718 // const MFloat dt = dgSolver().m_dt;
719 // initRestart(time, dt);
720 // initRestart(time, -1.0);
721 initRestart(time, m_fixedTimeStep); // NOTE: use prescribed fixed time step
722 } else {
723 // Load DG node variables and mean vars for source term calculation
724 calcInitialCondition(-1.0);
725 RECORD_TIMER_STOP(m_timers[Timers::Init]);
726
727 // Apply the initial conditions of the solver (again) to allow overwriting node variables set
728 // by the coupling
729 if(m_hasDgCartesianSolver) {
730 RECORD_TIMER_START(dgSolver().m_timers[maia::dg::Timers_::RunInit]);
731 RECORD_TIMER_START(dgSolver().m_timers[maia::dg::Timers_::InitData]);
732 dgSolver().initialCondition();
733 RECORD_TIMER_STOP(dgSolver().m_timers[maia::dg::Timers_::InitData]);
734 RECORD_TIMER_STOP(dgSolver().m_timers[maia::dg::Timers_::RunInit]);
735 }
736
737 RECORD_TIMER_START(m_timers[Timers::Init]);
738
739 // Update node variables at surfaces in the solver and perform extension
740 if(m_hasDgCartesianSolver && SysEqn::noNodeVars() > 0) {
741 dgSolver().updateNodeVariables();
742 dgSolver().extendNodeVariables();
743 }
744 }
745
746 RECORD_TIMER_STOP(m_timers[Timers::Init]);
747}
748
757template <MInt nDim, class CouplingDonor>
758void CouplingDgApe<nDim, CouplingDonor>::initProjection() {
759 TRACE();
760
761 if(!m_hasDgCartesianSolver) {
762 // Clear member variables if DG solver is not active
763 m_projection.clear();
764 m_elementDonorMap.clear();
765 m_elementDonorLevels.clear();
766 m_elementDonorPos.clear();
767 return;
768 }
769
770 // Create map from donor-cells to dg-elements
771 std::vector<std::vector<MInt>> elementGridMap(noElements());
772
773 for(MInt elementId = 0; elementId < noElements(); elementId++) {
774 const MInt cellId = elements().cellId(elementId);
775 const MInt gridCellId = dgSolver().grid().tree().solver2grid(cellId);
776
777 elementGridMap[elementId].clear();
778 // Get all leaf-cells of the donor-solver that are mapped to the global cell of the current
779 // dg-element
780 dgSolver().grid().raw().createLeafCellMapping(donorSolver().solverId(), gridCellId, elementGridMap[elementId]);
781 }
782
783 // Initialize the required Galerkin projections needed for spatial coupling
784 {
785 m_projection.clear();
786 m_projection.reserve(maxPolyDeg() + 1);
787
788 // Determine the number of refinement levels of the donor grid
789 MInt noLevels = 1;
790 if(m_hasDonorCartesianSolver) {
791 noLevels = donorSolver().grid().maxLevel() - donorSolver().grid().minLevel() + 1;
792 }
793
794 m_log << "Initialize Galerkin projections for " << noLevels << " levels." << std::endl;
795
796 // Dummy projections such that m_projection can be accessed by polyDeg
797 for(MInt i = 0; i < std::max(0, minPolyDeg() - m_noCutModesLowPass); i++) {
798 m_projection.push_back(ProjectionType());
799 }
800
801 // Initialize projection for each used polynomial degree
802 for(MInt i = std::max(0, minPolyDeg() - m_noCutModesLowPass); i <= maxPolyDeg(); i++) {
803 m_log << "Initialize Galerkin projection (p = " << i << ")... ";
804 m_projection.push_back(ProjectionType(i, noLevels, dgSolver().grid().lengthLevel0(), dgSolver().solverId()));
805 m_log << "done" << std::endl;
806 }
807 }
808
809 m_elementDonorMap.resize(noElements());
810 m_elementDonorLevels.resize(noElements());
811 m_elementDonorPos.resize(noElements());
812
813 // Determine information needed to apply the Galerkin projection
814 for(MInt elementId = 0; elementId < noElements(); elementId++) {
815 const MInt polyDeg = elements().polyDeg(elementId);
816
817 // Clear any previous projection information for this element
818 m_elementDonorMap[elementId].clear();
819 m_elementDonorLevels[elementId].clear();
820 m_elementDonorPos[elementId].clear();
821
822 // Check if element has mapped donor cells and skip iteration otherwise
823 if(elementGridMap[elementId].empty()) {
824 continue;
825 }
826
827 const MInt noMappedLeafCells = elementGridMap[elementId].size();
828 MIntScratchSpace donorCellLevels(noMappedLeafCells, AT_, "donorCellLevels");
829 MFloatScratchSpace donorCellCoordinates(noMappedLeafCells, nDim, AT_, "donorCellCoordinates");
830
831 // Information on target element
832 const MInt targetCellId = elements().cellId(elementId);
833 const MInt targetLevel = dgSolver().grid().tree().level(targetCellId);
834 const MFloat halfCellLength = 0.5 * dgSolver().grid().cellLengthAtLevel(targetLevel);
835 const MFloat* targetElementCenter = &dgSolver().grid().tree().coordinate(targetCellId, 0);
836
837 // Check if element is inside projection filter box (excluding filter)
838 if(m_projectionFilter) {
839 MBool hasProjection = false;
840 for(MInt dimId = 0; dimId < nDim; dimId++) {
841 const MFloat coord = targetElementCenter[dimId];
842 if(coord + halfCellLength < m_projectionFilterBox[dimId]
843 || coord - halfCellLength > m_projectionFilterBox[nDim + dimId]) {
844 hasProjection = true;
845 break;
846 }
847 }
848 if(!hasProjection) {
849 continue;
850 }
851 }
852
853 MInt noDonorCells = 0;
854 // Loop over all mapped donor leaf-cells
855 for(MInt i = 0; i < noMappedLeafCells; i++) {
856 const MInt donorGridCellId = elementGridMap[elementId][i];
857 const MInt donorCellId = donorSolver().grid().tree().grid2solver(donorGridCellId);
858
859 if(donorCellId < 0) {
860 TERMM(1, "Error in mapped donor leaf cells.");
861 }
862
863 // Store donor cell id in element donor map
864 m_elementDonorMap[elementId].push_back(donorCellId);
865 // Store the corresponding refinement level of the donor cell
866 donorCellLevels[noDonorCells] = donorSolver().grid().tree().level(donorCellId);
867
868 // Get original cell center coordinates
869 // NOTE: For all donor cells (i.e. also for boundary and cut cells) the
870 // original cell position is needed to establish the projection
871 // information. There is no special treatment of cut cells at the
872 // moment, for now they are treated as if they were normal cells.
873 for(MInt d = 0; d < nDim; d++) {
874 donorCellCoordinates(noDonorCells, d) = donorSolver().grid().tree().coordinate(donorCellId, d);
875 }
876
877 noDonorCells++;
878 }
879
880 if(noDonorCells == 0) {
881 TERMM(1, "Donor cells for element (elementId=" + std::to_string(elementId)
882 + ") not found, but elementGridMap[elementId] "
883 "with elementGridMap[elementId].size()="
884 + std::to_string(elementGridMap[elementId].size()) + " contains mapped donor cells.");
885 }
886
887 // Allocate storage for the projection information of this element
888 m_elementDonorLevels[elementId].resize(noDonorCells);
889 m_elementDonorPos[elementId].resize(noDonorCells);
890
891 // Check for one-to-multiple mapping and continue if so (there is no projection information to
892 // compute)
893 if(noDonorCells == 1) {
894 const MInt donorCellId = m_elementDonorMap[elementId][0];
895 const MInt donorLevel = donorSolver().grid().tree().level(donorCellId);
896 if(targetLevel > donorLevel) {
897 m_elementDonorLevels[elementId][0] = donorLevel - targetLevel;
898 m_elementDonorPos[elementId][0] = -1;
899 continue;
900 }
901 }
902
903 // Compute and store the projection information for this element (this includes non-mapped
904 // volumes)
905 m_projection[polyDeg].calcProjectionInformation(noDonorCells, &donorCellLevels[0], &donorCellCoordinates[0],
906 targetLevel, targetElementCenter, m_elementDonorPos[elementId],
907 m_elementDonorLevels[elementId]);
908 }
909}
910
917template <MInt nDim, class CouplingDonor>
918void CouplingDgApe<nDim, CouplingDonor>::initRestart(const MFloat time, const MFloat dt) {
919 TRACE();
920
921 // Store that this is a restart
922 m_isRestart = true;
923
924 // Calculate previous time for use in Qe source term
925 m_previousTime = time - dt;
926
927 RECORD_TIMER_START(m_timers[Timers::InitialCondition]);
928 // Only load mean data for source term calculation since node variables are already read from the
929 // nodevars-file
930 loadMeanQuantities(true);
931 RECORD_TIMER_STOP(m_timers[Timers::InitialCondition]);
932}
933
939template <MInt nDim, class CouplingDonor>
940void CouplingDgApe<nDim, CouplingDonor>::initSourceFilter() {
941 TRACE();
942
943 // Make sure that the element list is empty
944 m_calcSourceElements.clear();
945
946 if(!m_hasDgCartesianSolver) {
947 m_filter.clear();
948 m_elementInsideBody.clear();
949 return;
950 }
951
952 const MInt filterOffset = m_maxNoNodesXD;
953
954 // Temporary storage before permanent allocation with exact size
955 MBoolScratchSpace hasNonZeroFilter(noElements(), AT_, "hasNonZeroFilter");
956 std::fill(hasNonZeroFilter.begin(), hasNonZeroFilter.end(), false);
957
958 // Reset inside body marker for all elements
959 m_elementInsideBody.resize(noElements());
960 std::fill(m_elementInsideBody.begin(), m_elementInsideBody.end(), 0);
961
962 MInt noInsideElements = 0;
963
964 for(MInt elementId = 0; elementId < noElements(); elementId++) {
965 const MInt noNodes1D = elements().polyDeg(elementId) + 1;
966 const MInt noNodesXD = ipow(noNodes1D, nDim);
967
968 // Pointer to coordinates of source term
969 const MFloat* const coordinates = &elements().nodeCoordinates(elementId);
970
971 // Pointer to filter values of current element
972 MFloat* const elementFilter = &m_filter[elementId * filterOffset];
973
974 // Calculate filter values for all node coordinates of the element
975 m_sourceTermFilter.filter(coordinates, noNodesXD, elementFilter);
976
977 // Set filter values to zero if element is excluded via projection filter box
978 // @ansgar_mb TODO labels:COUPLER,DG make this a function?
979 if(m_projectionFilter) {
980 MBool hasProjection = false;
981 for(MInt dimId = 0; dimId < nDim; dimId++) {
982 const MFloat coord = coordinates[dimId];
983 if(coord < m_projectionFilterBox[dimId] || coord > m_projectionFilterBox[nDim + dimId]) {
984 hasProjection = true;
985 break;
986 }
987 }
988 if(!hasProjection) {
989 std::fill_n(elementFilter, noNodesXD, 0.0);
990 }
991 }
992
993 // Check if any filter value is nonzero (above epsilon), i.e. source terms need to be
994 // computed for this element
995 hasNonZeroFilter[elementId] =
996 std::any_of(elementFilter, elementFilter + noNodesXD, [](MFloat filter) { return filter > MFloatEps; });
997
998 // @ansgar TODO labels:COUPLER fix this, what about geometries that are not inside the source region -> mean
999 // vars (see loadMeanQuantities)?
1000 if(hasNonZeroFilter[elementId] && m_elementDonorMap[elementId].empty()) {
1001 m_elementInsideBody[elementId] = 1;
1002
1003 // Reset filter
1004 hasNonZeroFilter[elementId] = false;
1005 std::fill_n(elementFilter, noNodesXD, 0.0);
1006
1007 noInsideElements++;
1008 std::cerr << "Element " << std::to_string(elementId)
1009 << " has nonzero filter value(s), but there are no donor "
1010 "cells mapped to this element. It will be considered to be inside the geometry."
1011 << std::endl;
1012 TERMM(1, "TODO check! (element has nonzero filter but no donor cells)");
1013 }
1014 }
1015
1016 MPI_Allreduce(MPI_IN_PLACE, &noInsideElements, 1, maia::type_traits<MInt>::mpiType(), MPI_SUM, dgSolver().mpiComm(),
1017 AT_, "MPI_IN_PLACE", "noInsideElements");
1018 if(noInsideElements > 0) {
1019 m_log << "WARNING: Number of DG elements considered to be inside the geometry: " << noInsideElements << std::endl;
1020 if(dgSolver().domainId() == 0) {
1021 std::cerr << "WARNING: Number of DG elements considered to be inside the geometry: " << noInsideElements
1022 << std::endl;
1023 }
1024 }
1025
1026 // Permanently store elements with non-zero filter
1027 const MInt noNonZeroElements = std::count(hasNonZeroFilter.begin(), hasNonZeroFilter.end(), true);
1028 m_calcSourceElements.resize(noNonZeroElements);
1029
1030 // Store element id if source terms need to be computed, i.e. there is at
1031 // least one nonzero filter value AND there are donor cells mapped to this
1032 // element.
1033 for(MInt elementId = 0, nonZeroElementId = 0; elementId < noElements(); elementId++) {
1034 if(hasNonZeroFilter[elementId]) {
1035 if(!m_elementDonorMap[elementId].empty()) {
1036 m_calcSourceElements[nonZeroElementId] = elementId;
1037 nonZeroElementId++;
1038 } else {
1039 TERMM(1, "Element " + std::to_string(elementId)
1040 + " has nonzero filter value(s), but there are no donor "
1041 "cells mapped to this element. Go fix the source term "
1042 "filter in your property file.");
1043 }
1044 }
1045 }
1046
1047 // Compute spatial filter values on donor cells
1048 if(m_hasDonorCartesianSolver) {
1049 for(MInt i = 0; i < m_noDonorCells; i++) {
1050 const MFloat* const coordinates = &donorSolver().a_coordinate(i, 0);
1051 MFloat* const cellFilter = &m_filterDonor[i];
1052
1053 m_sourceTermFilter.filter(coordinates, 1, cellFilter);
1054 }
1055 }
1056}
1057
1062template <MInt nDim, class CouplingDonor>
1063void CouplingDgApe<nDim, CouplingDonor>::initTimers() {
1064 TRACE();
1065
1066 // Create timer group & timer for solver, and start the timer
1067 NEW_TIMER_GROUP_NOCREATE(m_timers[Timers::TimerGroup],
1068 "CouplingDgApe (solverId = " + std::to_string(solverId()) + ")");
1069 NEW_TIMER_NOCREATE(m_timers[Timers::Class], "total object lifetime", m_timers[Timers::TimerGroup]);
1070 RECORD_TIMER_START(m_timers[Timers::Class]);
1071
1072 // Create & start constructor timer
1073 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Constructor], "Constructor", m_timers[Timers::Class]);
1074 RECORD_TIMER_START(m_timers[Timers::Constructor]);
1075
1076 // Create regular solver-wide timers
1077 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Init], "init", m_timers[Timers::Class]);
1078 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::InitialCondition], "calcInitialCondition", m_timers[Timers::Class]);
1079 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::LoadMeanQuantities], "loadMeanQuantities",
1080 m_timers[Timers::InitialCondition]);
1081
1082 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::PreCouple], "PreCouple", m_timers[Timers::Class]);
1083 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::StoreSources], "storeSources", m_timers[Timers::PreCouple]);
1084 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::LoadInstantaneous], "loadInstantaneousQuantities",
1085 m_timers[Timers::StoreSources]);
1086
1087 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcSources], "calcSources", m_timers[Timers::StoreSources]);
1088 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcSourceLamb], "calcSourceLamb", m_timers[Timers::CalcSources]);
1089 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcSourceQmII], "calcSourceQmII", m_timers[Timers::CalcSources]);
1090 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcSourceQmIII], "calcSourceQmIII", m_timers[Timers::CalcSources]);
1091 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcSourceQe], "calcSourceQe", m_timers[Timers::CalcSources]);
1092 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::CalcSourceQc], "calcSourceQc", m_timers[Timers::CalcSources]);
1093
1094 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::ProjectSourceTerms], "projectSourceTerms", m_timers[Timers::StoreSources]);
1095 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::ApplyStoredSources], "applyStoredSources", m_timers[Timers::PreCouple]);
1096
1097 // Create accumulated timers that monitor selected subsystems
1098 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::Accumulated], "selected accumulated timers", m_timers[Timers::Class]);
1099 NEW_SUB_TIMER_NOCREATE(m_timers[Timers::LoadCouplingData], "loadCouplingData", m_timers[Timers::Accumulated]);
1100}
1101
1112template <MInt nDim, class CouplingDonor>
1113MBool CouplingDgApe<nDim, CouplingDonor>::loadCouplingData(const MString& filename,
1114 const MString& name_,
1115 const MInt stride,
1116 MFloat* data) {
1117 TRACE();
1118 RECORD_TIMER_START(m_timers[Timers::Accumulated]);
1119 RECORD_TIMER_START(m_timers[Timers::LoadCouplingData]);
1120
1121 using namespace maia::parallel_io;
1122 ParallelIo file(filename, PIO_READ, dgSolver().grid().raw().mpiComm());
1123
1124 MString datasetName;
1125 MString variableName;
1126 MBool foundVariable = false;
1127 // Get the names of all datasets in the file
1128 std::vector<MString> datasets = file.getDatasetNames();
1129 // Loop over all datasets and check if it has the right 'name' attribute
1130 for(auto&& dataset : datasets) {
1131 if(file.hasAttribute("name", dataset)) {
1132 file.getAttribute(&variableName, "name", dataset);
1133 if(variableName == name_) {
1134 datasetName = dataset;
1135 foundVariable = true;
1136 break;
1137 }
1138 }
1139 }
1140 // Read if variable was found in file
1141 if(foundVariable) {
1142 ParallelIo::size_type count = 0;
1143 ParallelIo::size_type start = 0;
1144
1145 if(m_hasDonorCartesianSolver) {
1146 count = donorSolver().noInternalCells();
1147 start = donorSolver().domainOffset(donorSolver().domainId());
1148 }
1149
1150 file.setOffset(count, start);
1151 file.readArray(data, datasetName, stride);
1152
1153 const MInt noInternalCells = donorSolver().noInternalCells();
1154 performUnitConversion(variableName, noInternalCells, stride, data);
1155 }
1156
1157 RECORD_TIMER_STOP(m_timers[Timers::LoadCouplingData]);
1158 RECORD_TIMER_STOP(m_timers[Timers::Accumulated]);
1159 return foundVariable;
1160}
1161
1166template <MInt nDim, class CouplingDonor>
1167void CouplingDgApe<nDim, CouplingDonor>::loadMeanQuantities(const MBool skipNodeVars) {
1168 TRACE();
1169 RECORD_TIMER_START(m_timers[Timers::LoadMeanQuantities]);
1170
1171 // Get the default node variables for regions without LES
1172 std::array<MFloat, SysEqn::noNodeVars()> defaultNodeVars;
1173 sysEqn().getDefaultNodeVars(&defaultNodeVars[0]);
1174
1175 // Get the default node variables for regions inside a LES geometry
1176 std::array<MFloat, SysEqn::noNodeVars()> defaultNodeVarsBody;
1177 sysEqn().getDefaultNodeVarsBody(&defaultNodeVarsBody[0]);
1178
1179 // Skip loading of the mean file (e.g. when doing a scaling and there is no useful mean file to be loaded)
1180 MBool skipLoadMeanFile = false;
1181 skipLoadMeanFile = Context::getSolverProperty<MBool>("skipLoadMeanFile", solverId(), AT_, &skipLoadMeanFile);
1182
1183 // Load mean data for node variables, apply the projection and store results
1184 // in nodeVars
1185 if(!skipNodeVars && !skipLoadMeanFile) {
1186 m_log << "Coupling condition: loading node vars from mean file." << std::endl;
1187 MFloatScratchSpace data(std::max(m_noDonorCells, 1), SysEqn::noNodeVars(), AT_, "data");
1188 // Load mean velocity data (data is stored with stride SysEqn::noNodeVars())
1189 std::vector<MInt> undefinedNodeVarIds; // nodeVarNames that are not available in mean file
1190 for(MInt varId = 0; varId < SysEqn::noNodeVars(); varId++) {
1191 if(!loadCouplingData(m_meanDataFileName, MV::nodeVarNames[varId], SysEqn::noNodeVars(), &data(0, varId))) {
1192 // Variable name was not found in mean file
1193 undefinedNodeVarIds.push_back(varId);
1194 }
1195 }
1196
1197 // Sanity check whether all data is available
1198 if(!undefinedNodeVarIds.empty()) {
1199 std::stringstream ss;
1200 ss << "Warning: Following dataset/s with attribute 'name' have not been found in file '" << m_meanDataFileName
1201 << "': ";
1202 for(MInt i : undefinedNodeVarIds) {
1203 ss << MV::nodeVarNames[i] << ", ";
1204 }
1205 ss << std::endl << "Those dataset/s have to be specified by the initial condition.";
1206 // Note: Here is only a warning called, since it might make sense not
1207 // providing some of the data through the mean file. For example the speed
1208 // of sound is not provided by the LB solver but set by the initial
1209 // condition as constant.
1210 // TERMM(1, ss.str());
1211 m_log << ss.str() << std::endl;
1212 }
1213
1214 // Disable low-pass filter for node vars
1215 const MBool previousLowPassFilterState = m_useLowPassFilter;
1216 m_useLowPassFilter = false;
1217
1218#ifdef _OPENMP
1219#pragma omp parallel for
1220#endif
1221 // Project data to elements and store in nodeVars
1222 for(MInt elementId = 0; elementId < noElements(); elementId++) {
1223 const MInt noDonorCells = m_elementDonorMap[elementId].size();
1224 const MInt noNonMappedCells = m_elementDonorPos[elementId].size() - noDonorCells;
1225
1226 // FIXME labels:COUPLER @ansgar this is not correct for a general case where there are bodies
1227 // in the donor-solver that are located outside the source region as DG-elements inside these
1228 // objects are not identified to be inside a body!
1229 if((m_elementInsideBody[elementId] == 0) && noNonMappedCells == 0) {
1230 projectToElement(elementId, SysEqn::noNodeVars(), &data[0], &defaultNodeVars[0],
1231 &elements().nodeVars(elementId));
1232 } else {
1233 projectToElement(elementId, SysEqn::noNodeVars(), &data[0], &defaultNodeVarsBody[0],
1234 &elements().nodeVars(elementId));
1235 }
1236 }
1237
1238 // Enable low-pass filter again after interpolation of node vars
1239 m_useLowPassFilter = previousLowPassFilterState;
1240 } else {
1241 m_log << "Coupling condition: skipped loading node vars from mean file." << std::endl;
1242 }
1243
1244 // Load mean data for source term calculation on donor grid
1245 if(!skipLoadMeanFile) {
1246 MFloatScratchSpace meanVars(std::max(m_noDonorCells, 1), noMeanVars(), AT_, "meanVars");
1247 m_log << "Coupling condition noMeanVars: " << noMeanVars() << std::endl;
1248 // Load mean variables data (stored with stride noMeanVars() in meanVars)
1249 std::vector<MInt> undefinedMeanVarIds; // meanVarNames that are not available in mean file
1250 for(MInt varId = 0; varId < noMeanVars(); varId++) {
1251 m_log << "Load coupling data: " << varId << " " << MV::names[m_activeMeanVars[varId]] << std::endl;
1252 if(!loadCouplingData(m_meanDataFileName, MV::names[m_activeMeanVars[varId]], noMeanVars(), &meanVars[varId])) {
1253 // Variable name was not found in mean file
1254 undefinedMeanVarIds.push_back(varId);
1255 }
1256 }
1257
1258 // Sanity check whether all data is available
1259 if(!undefinedMeanVarIds.empty()) {
1260 std::stringstream ss;
1261 ss << "Warning: Following dataset/s with attribute 'name' have not been found in file '" << m_meanDataFileName
1262 << "': ";
1263 for(MInt i : undefinedMeanVarIds) {
1264 ss << MV::names[m_activeMeanVars[i]] << ", ";
1265 }
1266 ss << std::endl;
1267 TERMM(1, ss.str());
1268 }
1269
1270 // Store the mean variables data of all donor leaf cells that are mapped to
1271 // elements that need to calculate coupling source terms in m_meanVars
1272 for(MInt donorIndex = 0; donorIndex < m_noActiveDonorCells; donorIndex++) {
1273 const MInt donorId = m_calcSourceDonorCells[donorIndex];
1274 std::copy_n(&meanVars(donorId, 0), noMeanVars(), &m_meanVars[donorIndex * noMeanVars()]);
1275 }
1276 } else {
1277 m_log << "Coupling condition: skip loading mean vars, set to zero instead." << std::endl;
1278 std::fill(m_meanVars.begin(), m_meanVars.end(), 0.0);
1279 }
1280
1281 RECORD_TIMER_STOP(m_timers[Timers::LoadMeanQuantities]);
1282}
1283
1296template <MInt nDim, class CouplingDonor>
1297void CouplingDgApe<nDim, CouplingDonor>::projectToElement(const MInt elementId, const MInt noVars, const MFloat* data,
1298 const MFloat* defaultValues, MFloat* target) {
1299 TRACE();
1300 const MInt polyDeg = elements().polyDeg(elementId);
1301 const MInt noNodesXD = ipow(polyDeg + 1, nDim);
1302 const MInt noDonorCells = m_elementDonorMap[elementId].size();
1303 // Number of non-mapped cells/volumes
1304 const MInt noNonMappedCells = m_elementDonorPos[elementId].size() - noDonorCells;
1305
1306 // Check if more than one donor cell
1307 if(noDonorCells > 1) {
1308 // Apply the Galerkin projection
1309 if(!m_useLowPassFilter || m_noCutModesLowPass == 0) {
1310 m_projection[polyDeg].apply(noDonorCells, noNonMappedCells, &m_elementDonorMap[elementId][0],
1311 &m_elementDonorLevels[elementId][0], &m_elementDonorPos[elementId][0], &data[0],
1312 defaultValues, noVars, target);
1313 } else {
1314 MFloatScratchSpace tmp(noNodesXD * noVars, AT_, "tmp");
1315 m_projection[polyDeg - m_noCutModesLowPass].apply(
1316 noDonorCells, noNonMappedCells, &m_elementDonorMap[elementId][0], &m_elementDonorLevels[elementId][0],
1317 &m_elementDonorPos[elementId][0], &data[0], defaultValues, noVars, &tmp[0]);
1318 maia::dg::interpolation::interpolateNodes<nDim>(&tmp[0],
1319 &m_vdmLowPass[polyDeg - m_noCutModesLowPass][0],
1320 polyDeg + 1 - m_noCutModesLowPass,
1321 polyDeg + 1,
1322 noVars,
1323 target);
1324 }
1325
1326 if(m_checkConservation || m_calcProjectionError) {
1327 if(noNonMappedCells > 0) {
1328 TERMM(1, "Fix conservation/projection error comp. for the case with non-mapped cells");
1329 }
1330
1331 const MInt cellId = elements().cellId(elementId);
1332 const MFloat targetLength = dgSolver().grid().cellLengthAtCell(cellId);
1333 MFloatScratchSpace errors(noVars, AT_, "errors");
1334
1335 // Check if the projection is conservative
1336 if(m_checkConservation) {
1337 m_projection[polyDeg].calcConservationError(noDonorCells, &m_elementDonorMap[elementId][0],
1338 &m_elementDonorLevels[elementId][0], &data[0], noVars, targetLength,
1339 target, &errors[0]);
1340
1341 const MFloat maxError = *(std::max_element(errors.begin(), errors.end()));
1342 m_maxConservationError = std::max(m_maxConservationError, maxError);
1343 }
1344
1345 // Calculate the projection error
1346 if(m_calcProjectionError) {
1347 m_projection[polyDeg].calcL2Error(noDonorCells, &m_elementDonorMap[elementId][0],
1348 &m_elementDonorLevels[elementId][0], &m_elementDonorPos[elementId][0],
1349 &data[0], noVars, targetLength, target, &errors[0]);
1350
1351 // Calculate maximum error
1352 const MFloat maxError = *(std::max_element(errors.begin(), errors.end()));
1353 m_maxL2Error = std::max(m_maxL2Error, maxError);
1354 }
1355 }
1356 } else if(noDonorCells == 1) {
1357 // One-to-one mapping
1358 MFloatTensor targetF(target, noNodesXD, noVars);
1359
1360 // Store donor cell variables in all element nodes
1361 const MInt donorCellId = m_elementDonorMap[elementId][0];
1362 const MInt offset = donorCellId * noVars;
1363 for(MInt i = 0; i < noNodesXD; i++) {
1364 for(MInt v = 0; v < noVars; v++) {
1365 targetF(i, v) = data[offset + v];
1366 }
1367 }
1368 } else {
1369 // No mapped cell: use default value(s)
1370 MFloatTensor targetF(target, noNodesXD, noVars);
1371 for(MInt i = 0; i < noNodesXD; i++) {
1372 for(MInt v = 0; v < noVars; v++) {
1373 targetF(i, v) = defaultValues[v];
1374 }
1375 }
1376 }
1377}
1378
1383template <MInt nDim, class CouplingDonor>
1384void CouplingDgApe<nDim, CouplingDonor>::saveFilterValues() {
1385 TRACE();
1386
1387 m_log << "Saving file with filter values ... ";
1388
1389 if(!m_hasDgCartesianSolver) {
1390 return;
1391 }
1392
1393 const std::vector<MString> filterVarName{"sourcefilter:" + m_sourceTermFilter.name()};
1394 saveNodalData("sourcefilter", 1, filterVarName, &m_filter[0]);
1395
1396 m_log << "done" << std::endl;
1397}
1398
1400template <MInt nDim, class CouplingDonor>
1401void CouplingDgApe<nDim, CouplingDonor>::saveFilterValuesDonor() {
1402 TRACE();
1403 m_log << "Saving file with filter values on donor cells... ";
1404
1405 if(!m_hasDonorCartesianSolver) {
1406 return;
1407 }
1408
1409 std::stringstream fileName;
1410 fileName << donorSolver().outputDir() << "sourcefilterDonor" << ParallelIo::fileExt();
1411
1412 using namespace maia::parallel_io;
1413 ParallelIo parallelIo(fileName.str(), PIO_REPLACE, donorSolver().mpiComm());
1414
1415 const MInt localNoCells = m_noDonorCells;
1416 // Determine offset and global number of cells
1417 ParallelIo::size_type offset = 0;
1418 ParallelIo::size_type globalNoCells = 0;
1419 parallelIo.calcOffset(localNoCells, &offset, &globalNoCells, donorSolver().mpiComm());
1420
1421 // Set attributes
1422 parallelIo.setAttribute(globalNoCells, "noCells");
1423 parallelIo.setAttribute(donorSolver().grid().gridInputFileName(), "gridFile");
1424 parallelIo.setAttribute(donorSolver().solverId(), "solverId");
1425
1426 const MString name_ = "sourcefilter:" + m_sourceTermFilter.name();
1427 parallelIo.defineArray(PIO_FLOAT, name_, globalNoCells);
1428 parallelIo.setOffset(localNoCells, offset);
1429 parallelIo.writeArray(&m_filterDonor[0], name_);
1430
1431 m_log << "done" << std::endl;
1432}
1433
1438template <MInt nDim, class CouplingDonor>
1439void CouplingDgApe<nDim, CouplingDonor>::applyStoredSources(const MFloat NotUsed(time)) {
1440 TRACE();
1441 RECORD_TIMER_START(m_timers[Timers::ApplyStoredSources]);
1442
1443 const MInt dataSize = m_maxNoNodesXD * SysEqn::noVars();
1444
1445 // Note: moved this into the DG solver applyExternalSource to remove the dependency
1446 // on the simulation time in this function, this way the external sources can be computed once in
1447 // the preCouple step before starting the first RK stage!
1448 // Source term ramp up factor (linear in time)
1449 /* const MFloat rampUpFactor */
1450 /* = (m_useSourceRampUp) */
1451 /* ? maia::filter::slope::linear(startTime(), */
1452 /* startTime() + m_sourceRampUpTime, time) */
1453 /* : 1.0; */
1454
1455 const MInt noSourceElements = m_calcSourceElements.size();
1456 // Add source terms from persistent storage to right-hand side of all
1457 // elements with both nonzero filter values and donor cells.
1458 maia::parallelFor(0, noSourceElements, [&](MInt elementIndex) {
1459 const MInt elementId = m_calcSourceElements[elementIndex];
1460 const MInt noNodes1D = elements().polyDeg(elementId) + 1;
1461 const MInt noNodes1D3 = (nDim == 3) ? noNodes1D : 1;
1462 MFloatTensor r(externalSource() + elementId * dataSize, noNodes1D, noNodes1D, noNodes1D3, SysEqn::noVars());
1463 MFloatTensor s(&m_sourceTerms[0] + elementIndex * dataSize, noNodes1D, noNodes1D, noNodes1D3, SysEqn::noVars());
1464
1465 // Add source terms to all equations
1466 for(MInt i = 0; i < noNodes1D; i++) {
1467 for(MInt j = 0; j < noNodes1D; j++) {
1468 for(MInt k = 0; k < noNodes1D3; k++) {
1469 for(MInt l = 0; l < SysEqn::noVars(); l++) {
1470 r(i, j, k, l) += m_sourceFactor * s(i, j, k, l);
1471 }
1472 }
1473 }
1474 }
1475 });
1476
1477 RECORD_TIMER_STOP(m_timers[Timers::ApplyStoredSources]);
1478}
1479
1491template <MInt nDim, class CouplingDonor>
1492void CouplingDgApe<nDim, CouplingDonor>::storeSources(const MFloat time, const MInt timeStep) {
1493 TRACE();
1494 RECORD_TIMER_START(m_timers[Timers::StoreSources]);
1495
1496 // Storage for velocities and vorticities of the donor cells
1497 MFloatScratchSpace velocity(std::max(m_noDonorCells, 1), noVelocities(), AT_, "velocity");
1498 MFloatScratchSpace vorticity(std::max(m_noDonorCells, 1), noVorticities(), AT_, "vorticity");
1499
1500 // Load all needed coupling data here
1501 // NOTE: for now (Q_mI and Q_mI_linear) the velocities and the vorticities are
1502 // needed for both source terms so there is no check yet which variables are
1503 // actually needed, but this will change in the future.
1504 RECORD_TIMER_START(m_timers[Timers::LoadInstantaneous]);
1505 getDonorVelocityAndVorticity(m_calcSourceDonorCells, velocity, vorticity);
1506 RECORD_TIMER_STOP(m_timers[Timers::LoadInstantaneous]);
1507
1508 // Source terms on donor cells (i.e. all source terms summed up)
1509 MFloatScratchSpace sourceTerms(std::max(m_noDonorCells, 1), SysEqn::noVars(), AT_, "sourceTerms");
1510 std::fill(sourceTerms.begin(), sourceTerms.end(), 0.0);
1511
1512 RECORD_TIMER_START(m_timers[Timers::CalcSources]);
1513 // Compute all source terms and accumulate to get the total source terms
1514 for(auto&& sourceTerm : m_activeSourceTerms) {
1515 switch(sourceTerm) {
1516 case ST::Q_mI:
1517 // Perturbed Lamb vector
1518 calcSourceLambNonlinear(&velocity[0], &vorticity[0], &sourceTerms(0, 0));
1519 break;
1520 case ST::Q_mI_linear:
1521 // Linearized Lamb vector
1522 calcSourceLambLinearized(&velocity[0], &vorticity[0], &sourceTerms(0, 0));
1523 break;
1524 case ST::Q_mII:
1525 // Pressure linearization
1526 calcSourceQmII(&velocity[0], &sourceTerms(0, 0));
1527 break;
1528 case ST::Q_mIII:
1529 calcSourceQmIII(&velocity[0], &sourceTerms(0, 0));
1530 break;
1531 case ST::Q_e:
1532 // Perturbed energy
1533 calcSourceQe(&velocity[0], time, &sourceTerms(0, 0));
1534 break;
1535 case ST::Q_c:
1536 calcSourceQc(&velocity[0], &sourceTerms(0, 0), time, timeStep);
1537 break;
1538 default:
1539 TERMM(1, "Source term '" + s_sourceTermNames[sourceTerm] + "' not implemented.");
1540 break;
1541 }
1542 }
1543 RECORD_TIMER_STOP(m_timers[Timers::CalcSources]);
1544
1545 // Set flag to indicate that the source terms have already been calculated at least once
1546 m_isFirstSourceCalculation = false;
1547
1548 // Store source terms on donor grid
1549 if(m_saveSourceTermsInterval > 0 && timeStep % m_saveSourceTermsInterval == 0 && m_saveSourceTermsDonorGrid) {
1550 stopLoadTimer(AT_);
1551 maia::dlb::g_dlbTimerController.disableAllDlbTimers();
1552
1553 saveSourceTermsDonorGrid(timeStep, &sourceTerms[0]);
1554
1555 maia::dlb::g_dlbTimerController.enableAllDlbTimers();
1556 startLoadTimer(AT_);
1557 }
1558
1559 // TODO labels:COUPLER save filtered or unfiltered sources on donor grid?
1560 if(m_applySourceFilterDonor) {
1561 // Apply spatial filter to source terms on donor cells
1562 const MInt noVars = SysEqn::noVars();
1563 for(MInt i = 0; i < m_noDonorCells; i++) {
1564 const MFloat filter = m_filterDonor[i];
1565 for(MInt j = 0; j < noVars; j++) {
1566 sourceTerms[i * noVars + j] *= filter;
1567 }
1568 }
1569 }
1570
1571 // Default source terms (NaN): only needed for the call to projectToElement()
1572 // but this should never be used/needed as an element without mapped LES cells
1573 // should not do anything here.
1574 /* const MFloat nan_value = std::numeric_limits<MFloat>::quiet_NaN(); */
1575 /* const std::array<MFloat, MAX_SPACE_DIMENSIONS + 1> defaultSourceTerms = { */
1576 /* {nan_value, nan_value, nan_value, nan_value}}; */
1577 // Note: changed to 0 since non-mapped cells can now appear in the projection information that
1578 // will use these values
1579 const std::array<MFloat, MAX_SPACE_DIMENSIONS + 1> defaultSourceTerms = {{0.0, 0.0, 0.0, 0.0}};
1580
1581 // Project accumulated source terms on elements and apply the spatial filter
1582 RECORD_TIMER_START(m_timers[Timers::ProjectSourceTerms]);
1583 const MInt noSourceElements = m_calcSourceElements.size();
1584 const MInt sourceOffset = m_maxNoNodesXD * SysEqn::noVars();
1585
1586 maia::parallelFor(0, noSourceElements, [&](MInt elementIndex) {
1587 const MInt elementId = m_calcSourceElements[elementIndex];
1588 const MInt noNodes1D = elements().polyDeg(elementId) + 1;
1589 const MInt noNodesXD = ipow(noNodes1D, nDim);
1590 const MInt elementOffset = elementIndex * sourceOffset;
1591
1592 MFloatTensor elementSourceTerms(&m_sourceTerms[elementOffset], noNodesXD, SysEqn::noVars());
1593 MFloatTensor filter(&m_filter[elementId * m_maxNoNodesXD], noNodesXD);
1594
1595 // Project source terms on element
1596 projectToElement(elementId, SysEqn::noVars(), &sourceTerms[0], &defaultSourceTerms[0], &elementSourceTerms[0]);
1597
1598 // Apply spatial filter on element (if not already applied on donor cells)
1599 if(!m_applySourceFilterDonor) {
1600 for(MInt i = 0; i < noNodesXD; i++) {
1601 for(MInt v = 0; v < SysEqn::noVars(); v++) {
1602 elementSourceTerms(i, v) *= filter(i);
1603 }
1604 }
1605 }
1606 });
1607 RECORD_TIMER_STOP(m_timers[Timers::ProjectSourceTerms]);
1608
1609 // Store previous time for calculating dp/dt in the Qe source term
1610 m_previousTime = time;
1611
1612 // Store source terms on target grid
1613 if(m_saveSourceTermsInterval > 0 && timeStep % m_saveSourceTermsInterval == 0) {
1614 stopLoadTimer(AT_);
1615 maia::dlb::g_dlbTimerController.disableAllDlbTimers();
1616
1617 saveSourceTermsTargetGrid(timeStep);
1618
1619 maia::dlb::g_dlbTimerController.enableAllDlbTimers();
1620 startLoadTimer(AT_);
1621 }
1622
1623 RECORD_TIMER_STOP(m_timers[Timers::StoreSources]);
1624}
1625
1633template <MInt nDim, class CouplingDonor>
1634void CouplingDgApe<nDim, CouplingDonor>::saveSourceTermsDonorGrid(const MInt timeStep, const MFloat* const data) {
1635 TRACE();
1636 CHECK_TIMERS_IO();
1637
1638 if(!m_hasDonorCartesianSolver) {
1639 return;
1640 }
1641
1642 const MInt noVariables = SysEqn::noVars();
1643
1644 std::stringstream fileName;
1645 fileName << donorSolver().outputDir() << "sourceTermsDonorGrid_" << std::setw(8) << std::setfill('0') << timeStep
1646 << ParallelIo::fileExt();
1647
1648 // Define source term names
1649 std::vector<MString> sourceTermNames(noVariables);
1650 if constexpr(nDim == 2) {
1651 sourceTermNames = std::vector<MString>{"source_u", "source_v", "source_p"};
1652 } else {
1653 sourceTermNames = std::vector<MString>{"source_u", "source_v", "source_w", "source_p"};
1654 }
1655
1656 // Create data out object to save data to disk
1657 using namespace maia::parallel_io;
1658 ParallelIo parallelIo(fileName.str(), PIO_REPLACE, donorSolver().mpiComm());
1659
1660 const MInt localNoCells = m_noDonorCells;
1661 // Determine offset and global number of cells
1662 ParallelIo::size_type offset = 0;
1663 ParallelIo::size_type globalNoCells = 0;
1664 parallelIo.calcOffset(localNoCells, &offset, &globalNoCells, donorSolver().mpiComm());
1665
1666 // Set attributes
1667 parallelIo.setAttribute(globalNoCells, "noCells");
1668 parallelIo.setAttribute(donorSolver().grid().gridInputFileName(), "gridFile");
1669 parallelIo.setAttribute(donorSolver().solverId(), "solverId");
1670
1671 // Define arrays in file
1672 for(MInt i = 0; i < noVariables; i++) {
1673 const MString name_ = "variables" + std::to_string(i);
1674 parallelIo.defineArray(PIO_FLOAT, name_, globalNoCells);
1675 parallelIo.setAttribute(sourceTermNames[i], "name", name_);
1676 }
1677
1678 parallelIo.setOffset(localNoCells, offset);
1679
1680 // Write source term data to file
1681 for(MInt i = 0; i < noVariables; i++) {
1682 const MString name_ = "variables" + std::to_string(i);
1683 parallelIo.writeArray(&data[i], name_, noVariables);
1684 }
1685}
1686
1693template <MInt nDim, class CouplingDonor>
1694void CouplingDgApe<nDim, CouplingDonor>::saveSourceTermsTargetGrid(const MInt timeStep) {
1695 TRACE();
1696 CHECK_TIMERS_IO();
1697
1698 const MInt noVariables = SysEqn::noVars();
1699 const MInt dataSize = m_maxNoNodesXD * noVariables;
1700
1701 // Buffer for source terms on all elements to write to file
1702 ScratchSpace<MFloat> sourceTerms(noElements(), dataSize, AT_, "sourceTerms");
1703
1704 // Copy source terms to buffer
1705 const MInt noSourceElements = m_calcSourceElements.size();
1706 for(MInt eId = 0; eId < noSourceElements; eId++) {
1707 const MInt elementId = m_calcSourceElements[eId];
1708 const MInt noNodesXD = ipow(elements().polyDeg(elementId) + 1, nDim);
1709
1710 // Set source and destination pointers
1711 MFloat* const dest = &sourceTerms[elementId * dataSize];
1712 const MFloat* const src = &m_sourceTerms[eId * dataSize];
1713
1714 // Copy all source terms of this element
1715 std::copy_n(src, noNodesXD * noVariables, dest);
1716 }
1717
1718 std::stringstream fileNameBase;
1719 fileNameBase << "sourceTerms_" << std::setw(8) << std::setfill('0') << timeStep;
1720
1721 // Define source term names
1722 std::vector<MString> sourceTermNames(noVariables);
1723 if constexpr(nDim == 2) {
1724 sourceTermNames = std::vector<MString>{"source_u", "source_v", "source_p"};
1725 } else {
1726 sourceTermNames = std::vector<MString>{"source_u", "source_v", "source_w", "source_p"};
1727 }
1728
1729 // Save source terms to file
1730 saveNodalData(fileNameBase.str(), noVariables, sourceTermNames, &sourceTerms[0]);
1731}
1732
1733
1739template <MInt nDim, class CouplingDonor>
1740void CouplingDgApe<nDim, CouplingDonor>::readProperties() {
1741 TRACE();
1742
1750 m_meanDataFileName = Context::getSolverProperty<MString>("meanDataFileName", solverId(), AT_);
1751
1760 // Get the number of source terms
1761 const MInt noSourceTerms = Context::propertyLength("sourceTerms", solverId());
1762
1763 // Ensure that at least one source term is loaded
1764 if(noSourceTerms == 0) {
1765 TERMM(1, "No source term was specified.");
1766 }
1767
1768 MBool q_mI_active = false;
1769 // Loop over all given source terms
1770 for(MInt i = 0; i < noSourceTerms; i++) {
1771 const MString sourceTermName = Context::getSolverProperty<MString>("sourceTerms", solverId(), AT_, i);
1772
1773 // Find the source term name in the list of all source terms
1774 auto sourceTermIt = std::find(s_sourceTermNames.begin(), s_sourceTermNames.end(), sourceTermName);
1775
1776 // Exit if source term name not found
1777 if(sourceTermIt == s_sourceTermNames.end()) {
1778 TERMM(1, "The given source term '" + sourceTermName + "' was not found.");
1779 }
1780
1781 // Determine source term index
1782 const MInt sourceTerm = std::distance(s_sourceTermNames.begin(), sourceTermIt);
1783
1784 // Check if this source term is already active
1785 if(std::find(m_activeSourceTerms.begin(), m_activeSourceTerms.end(), sourceTerm) != m_activeSourceTerms.end()) {
1786 TERMM(1, "Given source term '" + sourceTermName + "' already present. Check your property file.");
1787 }
1788
1789 // Make sure that only one q_mI source term is active
1790 if(sourceTerm == ST::Q_mI || sourceTerm == ST::Q_mI_linear) {
1791 if(q_mI_active) {
1792 TERMM(1, "Source terms 'q_mI' and 'q_mI_linear' cannot be active at the same time.");
1793 }
1794 q_mI_active = true;
1795 }
1796
1797 // Valid source term, store in list of active source terms
1798 m_activeSourceTerms.push_back(sourceTerm);
1799 }
1800 // Sort all active source terms by id
1801 std::sort(m_activeSourceTerms.begin(), m_activeSourceTerms.end());
1802
1803 m_log << "Direct-hybrid CFD-CAA multi-solver coupling condition" << std::endl;
1804 m_log << "Activated coupling source terms:" << std::endl;
1805 for(auto&& sourceTerm : m_activeSourceTerms) {
1806 m_log << s_sourceTermNames[sourceTerm] << std::endl;
1807 }
1808 m_log << "Number of active source terms: " << noSourceTerms << std::endl;
1809
1822 m_saveSourceTermsInterval = -1;
1823 m_saveSourceTermsInterval =
1824 Context::getSolverProperty<MInt>("saveSourceTermsInterval", solverId(), AT_, &m_saveSourceTermsInterval);
1825
1834 m_saveSourceTermsDonorGrid = false;
1835 m_saveSourceTermsDonorGrid =
1836 Context::getSolverProperty<MBool>("saveSourceTermsDonorGrid", solverId(), AT_, &m_saveSourceTermsDonorGrid);
1837
1838 // Determine unique list of needed mean variables for all active source terms
1839 // and store in m_activeMeanVars
1840 std::set<MInt> neededMeanVars;
1841 for(auto&& sourceTerm : m_activeSourceTerms) {
1842 std::vector<MInt> sourceTermMeanVars;
1843 // Get the needed mean variables for the current source term
1844 neededMeanVarsForSourceTerm(sourceTerm, sourceTermMeanVars);
1845 // Insert into set of needed mean variables
1846 neededMeanVars.insert(sourceTermMeanVars.begin(), sourceTermMeanVars.end());
1847 }
1848 m_activeMeanVars.assign(neededMeanVars.begin(), neededMeanVars.end());
1849
1850 // Store mapping from mean variable to corresponding index in m_activeMeanVars
1851 // Initialize with -1, since non-active mean variables should never be used
1852 std::fill_n(m_meanVarsIndex.begin(), s_totalNoMeanVars, -1);
1853 MInt meanVarPos = 0;
1854 // Loop over all active mean variables and store their position
1855 for(MInt meanVar : m_activeMeanVars) {
1856 m_meanVarsIndex[meanVar] = meanVarPos;
1857 meanVarPos++;
1858 }
1859
1867 m_fixedTimeStep = Context::getSolverProperty<MFloat>("fixedTimeStep", solverId(), AT_, &m_fixedTimeStep);
1868
1869 // Initialize source term filter
1870 m_sourceTermFilter.init();
1871
1880 m_applySourceFilterDonor = false;
1881 m_applySourceFilterDonor =
1882 Context::getSolverProperty<MBool>("applySourceFilterDonor", solverId(), AT_, &m_applySourceFilterDonor);
1883
1892 m_saveSourceTermFilter = false;
1893 m_saveSourceTermFilter =
1894 Context::getSolverProperty<MBool>("saveSourceTermFilter", solverId(), AT_, &m_saveSourceTermFilter);
1895
1904 m_checkConservation = false;
1905 m_checkConservation = Context::getSolverProperty<MBool>("checkConservation", solverId(), AT_, &m_checkConservation);
1906
1916 m_calcProjectionError = false;
1917 m_calcProjectionError =
1918 Context::getSolverProperty<MBool>("calcProjectionError", solverId(), AT_, &m_calcProjectionError);
1919
1927 m_sourceFactor = 1.0;
1928 m_sourceFactor = Context::getSolverProperty<MFloat>("sourceFactor", solverId(), AT_, &m_sourceFactor);
1929
1937 m_noCutModesLowPass = 0;
1938 m_noCutModesLowPass = Context::getSolverProperty<MInt>("noCutModesLowPass", solverId(), AT_, &m_noCutModesLowPass);
1939
1948 m_projectionFilter = false;
1949 m_projectionFilter = Context::getSolverProperty<MBool>("projectionFilter", solverId(), AT_, &m_projectionFilter);
1950
1951 if(m_projectionFilter) {
1952 m_log << "Using projection filter box" << std::endl;
1953 m_projectionFilterBox.resize(2 * nDim);
1954 for(MInt i = 0; i < 2 * nDim; i++) {
1962 m_projectionFilterBox[i] = Context::getSolverProperty<MFloat>("projectionFilterBox", solverId(), AT_, i);
1963 m_log << "projection filter box " << i << " " << m_projectionFilterBox[i] << std::endl;
1964 }
1965 }
1966}
1967
1968//--Methods (private)-----------------------------------------------------------
1969
1977template <MInt nDim, class CouplingDonor>
1978void CouplingDgApe<nDim, CouplingDonor>::neededMeanVarsForSourceTerm(const MInt sourceTerm,
1979 std::vector<MInt>& meanVars) const {
1980 TRACE();
1981
1982 meanVars.clear();
1983 switch(sourceTerm) {
1984 case ST::Q_mI:
1985 // Mean Lamb vector
1986 for(MInt i = 0; i < nDim; i++) {
1987 meanVars.push_back(MV::LAMB0[i]);
1988 }
1989 break;
1990 case ST::Q_mI_linear:
1991 // Mean velocities
1992 for(MInt i = 0; i < nDim; i++) {
1993 meanVars.push_back(MV::UU0[i]);
1994 }
1995 // Mean vorticities
1996 for(MInt i = 0; i < noVorticities(); i++) {
1997 meanVars.push_back(MV::VORT0[i]);
1998 }
1999 break;
2000 case ST::Q_mII:
2001 // Mean density and pressure
2002 meanVars.push_back(MV::RHO0);
2003 meanVars.push_back(MV::P0);
2004 // Mean gradient of rho
2005 for(MInt i = 0; i < nDim; i++) {
2006 meanVars.push_back(MV::DRHO[i]);
2007 }
2008 // Mean gradient of p
2009 for(MInt i = 0; i < nDim; i++) {
2010 meanVars.push_back(MV::DP[i]);
2011 }
2012 // Mean of (gradient of p divided by rho)
2013 for(MInt i = 0; i < nDim; i++) {
2014 meanVars.push_back(MV::GRADPRHO[i]);
2015 }
2016 break;
2017 case ST::Q_mIII:
2018 // Mean velocities
2019 for(MInt i = 0; i < nDim; i++) {
2020 meanVars.push_back(MV::UU0[i]);
2021 }
2022 // Mean velocity gradients
2023 for(MInt i = 0; i < nDim * nDim; i++) {
2024 meanVars.push_back(MV::GRADU[i]);
2025 }
2026 // Sum of products of velocity components with corresponding velocity gradients
2027 for(MInt i = 0; i < nDim; i++) {
2028 meanVars.push_back(MV::UGRADU[i]);
2029 }
2030 break;
2031 case ST::Q_e:
2032 meanVars.push_back(MV::RHO0);
2033 meanVars.push_back(MV::P0);
2034 meanVars.push_back(MV::C0);
2035 // Mean velocities
2036 for(MInt i = 0; i < nDim; i++) {
2037 meanVars.push_back(MV::UU0[i]);
2038 }
2039 // Mean gradient of c
2040 for(MInt i = 0; i < nDim; i++) {
2041 meanVars.push_back(MV::DC0[i]);
2042 }
2043 // Components of divergence of u
2044 for(MInt i = 0; i < nDim; i++) {
2045 meanVars.push_back(MV::DU[i]);
2046 }
2047 // Mean gradient of rho
2048 for(MInt i = 0; i < nDim; i++) {
2049 meanVars.push_back(MV::DRHO[i]);
2050 }
2051 // Mean gradient of p
2052 for(MInt i = 0; i < nDim; i++) {
2053 meanVars.push_back(MV::DP[i]);
2054 }
2055 break;
2056 case ST::Q_c:
2057 // Mean density, pressure, and speed of sound
2058 meanVars.push_back(MV::RHO0);
2059 meanVars.push_back(MV::P0);
2060 meanVars.push_back(MV::C0);
2061 // Mean velocities
2062 for(MInt i = 0; i < nDim; i++) {
2063 meanVars.push_back(MV::UU0[i]);
2064 }
2065 // Components of divergence of u
2066 for(MInt i = 0; i < nDim; i++) {
2067 meanVars.push_back(MV::DU[i]);
2068 }
2069 // Mean gradient of rho
2070 for(MInt i = 0; i < nDim; i++) {
2071 meanVars.push_back(MV::DRHO[i]);
2072 }
2073 // Mean gradient of rho*div(u)
2074 for(MInt i = 0; i < nDim; i++) {
2075 meanVars.push_back(MV::RHODIVU[i]);
2076 }
2077 // Mean gradient of u*grad(rho)
2078 for(MInt i = 0; i < nDim; i++) {
2079 meanVars.push_back(MV::UGRADRHO[i]);
2080 }
2081 break;
2082 default:
2083 TERMM(1, "Source term '" + s_sourceTermNames[sourceTerm] + "' not implemented yet.");
2084 break;
2085 }
2086}
2087
2088#endif // COUPLINGDGAPE_H_
Context::propertyLength
static MInt propertyLength(const MString &name, MInt solverId=m_noSolvers)
Returns the number of elements of a property.
Definition: context.cpp:538
Context::getSolverProperty
static T getSolverProperty(const MString name, const MInt solverId, const MString &nameOfCallingFunction, const T *default_value, MInt pos=0)
Definition: context.h:168
CouplingDgApe
Intermediate class for coupling DG solver with APE sysEqn.
Definition: couplingdgape.h:214
CouplingDgApe::m_elementDonorMap
std::vector< std::vector< MInt > > m_elementDonorMap
Mapping from donor cells to elements.
Definition: couplingdgape.h:381
CouplingDgApe::postCouple
void postCouple(MInt) override final
Definition: couplingdgape.h:258
CouplingDgApe::m_elementDonorPos
std::vector< std::vector< MInt > > m_elementDonorPos
Donor cell positions on the corresponding cell level.
Definition: couplingdgape.h:382
CouplingDgApe::s_sourceTermNames
static const std::array< MString, ST::totalNoSourceTerms > s_sourceTermNames
Hold indices for source terms.
Definition: couplingdgape.h:430
CouplingDgApe::projectToElement
void projectToElement(const MInt elementId, const MInt noVars, const MFloat *data, const MFloat *defaultValues, MFloat *target)
Project the given data fields onto a single element.
Definition: couplingdgape.h:1297
CouplingDgApe::m_sourceFactor
MFloat m_sourceFactor
Factor by which the source terms are multiplied.
Definition: couplingdgape.h:392
CouplingDgApe::m_activeMeanVars
std::vector< MInt > m_activeMeanVars
List of active mean variables for all active source terms.
Definition: couplingdgape.h:369
CouplingDgApe::init
void init() override
Definition: couplingdgape.h:249
CouplingDgApe::performUnitConversion
virtual void performUnitConversion(const MString &, const MInt, const MInt, MFloat *)
Definition: couplingdgape.h:331
CouplingDgApe::m_filter
std::vector< MFloat > m_filter
Local storage for source filter values.
Definition: couplingdgape.h:408
CouplingDgApe::loadCouplingData
MBool loadCouplingData(const MString &filename, const MString &name_, const MInt stride, MFloat *data)
Auxiliary method to load coupling data (i.e. anything coming from the LES) from a file.
Definition: couplingdgape.h:1113
CouplingDgApe::m_vdmLowPass
std::vector< MFloatTensor > m_vdmLowPass
Vandermonde matrix/matrices for the low-pass filter.
Definition: couplingdgape.h:412
CouplingDgApe::noVelocities
static constexpr MInt noVelocities()
Return number of velocity variables.
Definition: couplingdgape.h:306
CouplingDgApe::m_sourceTerms
std::vector< MFloat > m_sourceTerms
Local storage for source terms.
Definition: couplingdgape.h:389
CouplingDgApe::getCouplingTimings
void getCouplingTimings(std::vector< std::pair< MString, MFloat > > &timings, const MBool allTimings) override
Definition: couplingdgape.h:261
CouplingDgApe::m_calcSourceElements
std::vector< MInt > m_calcSourceElements
List of all element ids for which source terms need to be computed.
Definition: couplingdgape.h:391
CouplingDgApe::m_isRestart
MBool m_isRestart
Store whether this is a restart (in case special treatment is necessary)
Definition: couplingdgape.h:356
CouplingDgApe::m_checkConservation
MBool m_checkConservation
Check if each Galerkin projection is conservative.
Definition: couplingdgape.h:395
CouplingDgApe::calcSourceLambLinearized
virtual void calcSourceLambLinearized(const MFloat *const velocity, const MFloat *const vorticity, MFloat *sourceTerms)=0
CouplingDgApe::saveSourceTermsDonorGrid
void saveSourceTermsDonorGrid(const MInt timeStep, const MFloat *const data)
Store the source terms on the donor grid.
Definition: couplingdgape.h:1634
CouplingDgApe::m_projectionFilter
MBool m_projectionFilter
Use spatial projection filter (i.e. exclude elements from the projection)
Definition: couplingdgape.h:416
CouplingDgApe::initSourceFilter
void initSourceFilter()
Calculate source filter values for all elements and store element ids on which source terms need to b...
Definition: couplingdgape.h:940
CouplingDgApe::m_meanVarsIndex
std::array< MInt, s_totalNoMeanVars > m_meanVarsIndex
Definition: couplingdgape.h:426
CouplingDgApe::m_isFirstSourceCalculation
MBool m_isFirstSourceCalculation
Store whether this is the first calculation of the source terms.
Definition: couplingdgape.h:388
CouplingDgApe::noMeanVars
MInt noMeanVars() const
Return number of mean variables.
Definition: couplingdgape.h:308
CouplingDgApe::m_elementDonorLevels
std::vector< std::vector< MInt > > m_elementDonorLevels
Donor cell levels relative to element.
Definition: couplingdgape.h:380
CouplingDgApe::m_filterDonor
std::vector< MFloat > m_filterDonor
Local storage for source filter values on donor cells.
Definition: couplingdgape.h:409
CouplingDgApe::m_previousTimeStep
MInt m_previousTimeStep
Previous time step (to determine whether new sources need to be calculated)
Definition: couplingdgape.h:367
CouplingDgApe::saveSourceTermsTargetGrid
void saveSourceTermsTargetGrid(const MInt timeStep)
Store the source terms on the target grid, i.e. after interpolation.
Definition: couplingdgape.h:1694
CouplingDgApe::m_sourceTermFilter
Filter< nDim > m_sourceTermFilter
Auxiliary object that handles source filtering.
Definition: couplingdgape.h:407
CouplingDgApe::initData
void initData()
Initialize the data (initial conditions or for a restart) of the coupling condition.
Definition: couplingdgape.h:710
CouplingDgApe::calcSourceQmIII
virtual void calcSourceQmIII(const MFloat *const velocity, MFloat *sourceTerms)=0
CouplingDgApe::loadMeanQuantities
void loadMeanQuantities(const MBool skipNodeVars=false)
Load mean velocities from file and store them to the node variables.
Definition: couplingdgape.h:1167
CouplingDgApe::saveFilterValues
void saveFilterValues()
Save filter values to file.
Definition: couplingdgape.h:1384
CouplingDgApe::m_timers
std::array< MInt, Timers::_count > m_timers
Definition: couplingdgape.h:464
CouplingDgApe::calcSourceQe
virtual void calcSourceQe(const MFloat *const velocity, const MFloat time, MFloat *const sourceTerms)=0
CouplingDgApe::initProjection
void initProjection()
Initialize the projection information for spatial coupling.
Definition: couplingdgape.h:758
CouplingDgApe::saveFilterValuesDonor
void saveFilterValuesDonor()
Save the filter values on the donor cells.
Definition: couplingdgape.h:1401
CouplingDgApe::noVorticities
static constexpr MInt noVorticities()
Return number of vorticity variables.
Definition: couplingdgape.h:307
CouplingDgApe::m_calcProjectionError
MBool m_calcProjectionError
Calculate the L2 error of the Galerkin projection.
Definition: couplingdgape.h:374
CouplingDgApe::calcSourceQc
virtual void calcSourceQc(const MFloat *const velocity, MFloat *const sourceTerms, const MFloat time, const MInt timeStep)=0
CouplingDgApe::CouplingDgApe
CouplingDgApe(const MInt couplingId, DgCartesianSolverType *dg, DonorSolverType *ds)
Initialize timers and read properties in c'tor.
Definition: couplingdgape.h:500
CouplingDgApe::donorSolver
virtual DonorSolverType & donorSolver(const MInt solverId=0) const =0
CouplingDgApe::m_meanDataFileName
MString m_meanDataFileName
File name for mean quantities.
Definition: couplingdgape.h:370
CouplingDgApe::applyStoredSources
void applyStoredSources(const MFloat time)
Apply the stored source terms to the time derivative.
Definition: couplingdgape.h:1439
CouplingDgApe::m_maxL2Error
MFloat m_maxL2Error
Maximum computed L2 error of the Galerkin projection.
Definition: couplingdgape.h:397
CouplingDgApe::m_calcSourceDonorCells
std::vector< MInt > m_calcSourceDonorCells
List of all donor cell ids for which source terms need to be computed.
Definition: couplingdgape.h:390
CouplingDgApe::m_saveSourceTermFilter
MBool m_saveSourceTermFilter
Store whether filter values should be written to a file at initialization.
Definition: couplingdgape.h:414
CouplingDgApe::m_applySourceFilterDonor
MBool m_applySourceFilterDonor
Apply source filter on donor cells or on DG elements after projection.
Definition: couplingdgape.h:413
CouplingDgApe::initRestart
void initRestart(const MFloat time, const MFloat dt)
Perform initializations for a restart.
Definition: couplingdgape.h:918
CouplingDgApe::m_projection
std::vector< ProjectionType > m_projection
Galerkin projection (access by polynomial degree)
Definition: couplingdgape.h:385
CouplingDgApe::m_noDonorCells
MInt m_noDonorCells
Total number of donor cells on this domain.
Definition: couplingdgape.h:379
CouplingDgApe::m_noCutModesLowPass
MInt m_noCutModesLowPass
Number of modes to cut using the low-pass source term.
Definition: couplingdgape.h:410
CouplingDgApe::m_saveSourceTermsInterval
MInt m_saveSourceTermsInterval
Interval at which source term data should be saved to disk.
Definition: couplingdgape.h:373
CouplingDgApe::getDonorVelocityAndVorticity
virtual void getDonorVelocityAndVorticity(const std::vector< MInt > &donorCellIds, MFloatScratchSpace &p_velocity, MFloatScratchSpace &p_vorticity)=0
CouplingDgApe::checkProperties
void checkProperties() override final
Definition: couplingdgape.h:348
CouplingDgApe::m_previousTime
MFloat m_previousTime
Previous time for the calculation of time derivatives.
Definition: couplingdgape.h:365
CouplingDgApe::initCoupler
void initCoupler()
Initialize the coupling condition (data structures).
Definition: couplingdgape.h:566
CouplingDgApe::m_noActiveDonorCells
MInt m_noActiveDonorCells
Definition: couplingdgape.h:377
CouplingDgApe::m_useLowPassFilter
MBool m_useLowPassFilter
Switch low pass on and off to allow disabling e.g. for node vars.
Definition: couplingdgape.h:411
CouplingDgApe::getDomainDecompositionInformation
void getDomainDecompositionInformation(std::vector< std::pair< MString, MInt > > &domainInfo) override final
Definition: couplingdgape.h:284
CouplingDgApe::m_hasDonorCartesianSolver
MBool m_hasDonorCartesianSolver
Store whether this domain has Cartesian donor solver cells.
Definition: couplingdgape.h:358
CouplingDgApe::finalizeSubCoupleInit
void finalizeSubCoupleInit(MInt) override final
Definition: couplingdgape.h:254
CouplingDgApe::preCouple
void preCouple(MInt) final
Calculate source terms and add to the external source terms of the DG solver.
Definition: couplingdgape.h:519
CouplingDgApe::calcSourceQmII
virtual void calcSourceQmII(const MFloat *const velocity, MFloat *const sourceTerms)=0
CouplingDgApe::BaseDonor
CouplingDonor BaseDonor
Definition: couplingdgape.h:221
CouplingDgApe::m_fixedTimeStep
MFloat m_fixedTimeStep
Fixed time step to use.
Definition: couplingdgape.h:364
CouplingDgApe::m_projectionFilterBox
std::vector< MFloat > m_projectionFilterBox
Spatial projection filter box (excluding)
Definition: couplingdgape.h:417
CouplingDgApe::m_activeSourceTerms
std::vector< MInt > m_activeSourceTerms
List of active source terms.
Definition: couplingdgape.h:371
CouplingDgApe::DonorSolverType
typename BaseDonor::solverType DonorSolverType
Definition: couplingdgape.h:222
CouplingDgApe::~CouplingDgApe
virtual ~CouplingDgApe()=default
CouplingDgApe::subCouple
void subCouple(MInt, MInt, std::vector< MBool > &) override final
Definition: couplingdgape.h:257
CouplingDgApe::initTimers
void initTimers()
Create timers.
Definition: couplingdgape.h:1063
CouplingDgApe::neededMeanVarsForSourceTerm
void neededMeanVarsForSourceTerm(const MInt sourceTerm, std::vector< MInt > &meanVars) const
Return the needed mean variables for a given source term.
Definition: couplingdgape.h:1978
CouplingDgApe::m_hasDgCartesianSolver
MBool m_hasDgCartesianSolver
Store whether this domain has DG cells/elements.
Definition: couplingdgape.h:357
CouplingDgApe::m_elementInsideBody
std::vector< MInt > m_elementInsideBody
Marker if elements are inside a geometric object.
Definition: couplingdgape.h:403
CouplingDgApe::noCouplingTimers
MInt noCouplingTimers(const MBool allTimings) const override
Definition: couplingdgape.h:293
CouplingDgApe::m_meanVars
std::vector< MFloat > m_meanVars
Local storage for mean variables of the donor cells.
Definition: couplingdgape.h:400
CouplingDgApe::cleanUp
void cleanUp() override final
Definition: couplingdgape.h:259
CouplingDgApe::m_maxConservationError
MFloat m_maxConservationError
Maximum conservation error.
Definition: couplingdgape.h:396
CouplingDgApe::storeSources
void storeSources(const MFloat time, const MInt timeStep)
Load coupling data from LES, then calculate, accumulate, project and store coupling source terms.
Definition: couplingdgape.h:1492
CouplingDgApe::m_maxNoNodesXD
MInt m_maxNoNodesXD
Maximum number of nodes of an element (corresponding to maxPolyDeg)
Definition: couplingdgape.h:360
CouplingDgApe::m_saveSourceTermsDonorGrid
MBool m_saveSourceTermsDonorGrid
Store whether the sources on the donor grid should be saved as well.
Definition: couplingdgape.h:372
CouplingDgApe::calcSourceLambNonlinear
virtual void calcSourceLambNonlinear(const MFloat *const velocity, const MFloat *const vorticity, MFloat *const sourceTerms)=0
CouplingDgApe::calcInitialCondition
void calcInitialCondition(const MFloat time)
Apply initial conditions for this coupling condition.
Definition: couplingdgape.h:549
CouplingDgApe::readProperties
void readProperties() override final
Read properties and store in member variables.
CouplingDgApe::s_totalNoMeanVars
static const MInt s_totalNoMeanVars
Definition: couplingdgape.h:422
CouplingDgApe::DgCartesianSolverType
typename BaseDg::solverType DgCartesianSolverType
Definition: couplingdgape.h:227
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::externalSource
MFloat * externalSource() const
Return pointer to external source memory.
Definition: coupling.h:327
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::solverType
DgCartesianSolver< nDim, DgSysEqnAcousticPerturb< nDim > > solverType
Definition: coupling.h:292
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::noElements
MInt noElements() const
Return number of elements.
Definition: coupling.h:324
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::minPolyDeg
MInt minPolyDeg() const
Return the minimum polynomial degree.
Definition: coupling.h:333
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::sysEqn
DgSysEqnAcousticPerturb< nDim > & sysEqn()
Return reference to SysEqn object.
Definition: coupling.h:318
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::saveNodalData
void saveNodalData(const MString &fileNameBase, const MInt noVars, const std::vector< MString > &varNames, const MFloat *const data) const
Save nodal data to file.
Definition: coupling.h:342
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::elements
ElementCollector & elements()
Return reference to elements.
Definition: coupling.h:321
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::dgSolver
solverType & dgSolver() const
Return MPI communicator.
Definition: coupling.h:310
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::solverId
MInt solverId() const
Return solver id.
Definition: coupling.h:315
CouplingDg< nDim, DgSysEqnAcousticPerturb< nDim > >::maxPolyDeg
MInt maxPolyDeg() const
Return the maximum polynomial degree.
Definition: coupling.h:336
Coupling::returnIdleRecord
MFloat returnIdleRecord() const
Definition: coupling.h:159
Coupling::stopLoadTimer
void stopLoadTimer(const MString &name) const
Stop the load timer of the coupler.
Definition: coupling.h:152
Coupling::startLoadTimer
void startLoadTimer(const MString &name) const
Start the load timer of the coupler.
Definition: coupling.h:148
Coupling::returnLoadRecord
MFloat returnLoadRecord() const
Definition: coupling.h:158
Coupling::couplerId
MInt couplerId() const
Definition: coupling.h:67
DgGalerkinProjection
Performs the Galerkin projection for a given polynomial degree.
Definition: dgcartesiangalerkinprojection.h:32
DgInterpolation
Class stores precalculated values for interpolation & integration on the reference interval [-1,...
Definition: dgcartesianinterpolation.h:35
DlbTimerController::enableAllDlbTimers
void enableAllDlbTimers(const MBool *const wasEnabled=nullptr)
Enable all DLB timers (or those given by the array wasEnabled)
Definition: dlbtimer.h:265
DlbTimerController::disableAllDlbTimers
void disableAllDlbTimers(MBool *const wasEnabled=nullptr)
Disable all (enabled) DLB timers.
Definition: dlbtimer.h:300
Filter
Filter object for source terms.
Definition: filter.h:171
ScratchSpace
This class is a ScratchSpace.
Definition: scratch.h:758
ScratchSpace::end
iterator end()
Definition: scratch.h:281
ScratchSpace::begin
iterator begin()
Definition: scratch.h:273
string2enum
MInt string2enum(MString theString)
This global function translates strings in their corresponding enum values (integer values)....
Definition: enums.cpp:20
DgPolynomialType
DgPolynomialType
Definition: enums.h:310
DgIntegrationMethod
DgIntegrationMethod
Definition: enums.h:313
ipow
MInt ipow(MInt base, MInt exp)
Integer exponent function for non-negative exponents.
Definition: functions.h:317
m_log
InfoOutFile m_log
Definition: globalvariables.cpp:57
MInt
int32_t MInt
Definition: maiatypes.h:62
MString
std::basic_string< char > MString
Definition: maiatypes.h:55
MFloat
double MFloat
Definition: maiatypes.h:52
MBool
bool MBool
Definition: maiatypes.h:58
MPI_Allreduce
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
Definition: mpioverride.cpp:952
maia::dg::interpolation::calcPolynomialInterpolationMatrix
void calcPolynomialInterpolationMatrix(MInt noNodesIn, const T nodesIn, MInt noNodesOut, const U nodesOut, const V wBary, W vandermonde)
Calculate the polynomial interpolation matrix (Vandermonde) to interpolate from one set of nodes to a...
Definition: dgcartesianinterpolation.h:729
maia::dlb::g_dlbTimerController
DlbTimerController g_dlbTimerController
Definition: globalvariables.cpp:42
maia::parallel_io
Definition: parallelio.h:26
maia::parallelFor
void parallelFor(MInt begin, MInt end, UnaryFunction &&f)
Wrapper function for parallel for loop.
Definition: parallelfor.h:147
ParallelIo
PARALLELIO_DEFAULT_BACKEND ParallelIo
Definition: parallelio.h:292
MeanVariables< 2 >::names
static const std::array< MString, totalNoMeanVars > names
Definition: couplingdgape.h:111
MeanVariables< 2 >::nodeVarNames
static const std::array< MString, 6 > nodeVarNames
Definition: couplingdgape.h:112
MeanVariables< 3 >::nodeVarNames
static const std::array< MString, 8 > nodeVarNames
Definition: couplingdgape.h:206
MeanVariables< 3 >::names
static const std::array< MString, totalNoMeanVars > names
Definition: couplingdgape.h:205
maia::dg::Timers_::InitData
@ InitData
Definition: dgcartesiantimers.h:24
maia::dg::Timers_::RunInit
@ RunInit
Definition: dgcartesiantimers.h:20
maia::type_traits
Definition: typetraits.h:22