Nektar++
VariableConverter.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: VariableConverter.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Auxiliary functions to convert variables in
32// the compressible flow system
33//
34///////////////////////////////////////////////////////////////////////////////
35
36#include <iomanip>
37#include <iostream>
38
43
44namespace Nektar
45{
46
48 const LibUtilities::SessionReaderSharedPtr &pSession, const int spaceDim,
50 : m_session(pSession), m_spacedim(spaceDim)
51{
52 // Create equation of state object
53 std::string eosType;
54 m_session->LoadSolverInfo("EquationOfState", eosType, "IdealGas");
56
57 // Parameters for dynamic viscosity
58 m_session->LoadParameter("pInf", m_pInf, 101325);
59 m_session->LoadParameter("rhoInf", m_rhoInf, 1.225);
60 m_session->LoadParameter("GasConstant", m_gasConstant, 287.058);
61 m_session->LoadParameter("mu", m_mu, 1.78e-05);
63
64 // Parameters for sensor
65 m_session->LoadParameter("Skappa", m_Skappa, -1.0);
66 m_session->LoadParameter("Kappa", m_Kappa, 0.25);
67
69
70 // Shock sensor
71 m_session->LoadSolverInfo("ShockCaptureType", m_shockCaptureType, "Off");
72 if (m_shockCaptureType == "Physical")
73 {
74 // Artificial viscosity scaling constant
75 m_session->LoadParameter("mu0", m_mu0, 1.0);
76
79
80 // Check for Modal/Dilatation sensor
81 m_session->LoadSolverInfo("ShockSensorType", m_shockSensorType,
82 "Dilatation");
83
84 // Check for Ducros sensor
85 m_session->LoadSolverInfo("DucrosSensor", m_ducrosSensor, "Off");
86
87 if (m_ducrosSensor != "Off" || m_shockSensorType == "Dilatation")
88 {
89 m_flagCalcDivCurl = true;
90 }
91 }
92 // Load smoothing type.
93 // We only allocate the smoother if a mesh-graph was passed to the
94 // constructor. This is done to prevent the smoother from being allocated
95 // in cases when it won't be used. Otherwise, any class in the code that
96 // allocates a VariableConverter object try to will allocate a smoother,
97 // even if no mesh graph was passed to the constructor.
98 // TODO: the smoother should be separated from the VariableConverter class.
99 m_session->LoadSolverInfo("Smoothing", m_smoothing, "Off");
100 if (m_smoothing == "C0" && pGraph.get() != nullptr)
101 {
104 m_session, pGraph, m_session->GetVariable(0));
105 }
106
107 std::string viscosityType;
108 m_session->LoadSolverInfo("ViscosityType", viscosityType, "Constant");
109 if ("Variable" == viscosityType)
110 {
111 WARNINGL0(
112 m_session->DefinesParameter("Tref"),
113 "The Tref should be given in Kelvin for using the Sutherland's law "
114 "of dynamic viscosity. The default is 288.15. Note the mu or "
115 "Reynolds number should coorespond to this temperature.");
116 m_session->LoadParameter("Tref", m_Tref, 288.15);
117 m_TRatioSutherland = 110.0 / m_Tref;
118 }
119}
120
121/**
122 * @brief Compute the dynamic energy
123 * \f$ e = rho*V^2/2 \f$.
124 */
126 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
128{
129 size_t nPts = physfield[m_spacedim + 1].size();
130 Vmath::Zero(nPts, energy, 1);
131
132 // tmp = (rho * u_i)^2
133 for (size_t i = 0; i < m_spacedim; ++i)
134 {
135 Vmath::Vvtvp(nPts, physfield[i + 1], 1, physfield[i + 1], 1, energy, 1,
136 energy, 1);
137 }
138 // Divide by rho and multiply by 0.5 --> tmp = 0.5 * rho * u^2
139 Vmath::Vdiv(nPts, energy, 1, physfield[0], 1, energy, 1);
140 Vmath::Smul(nPts, 0.5, energy, 1, energy, 1);
141}
142
143/**
144 * @brief Compute the specific internal energy
145 * \f$ e = (E - rho*V^2/2)/rho \f$.
146 */
148 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
150{
151 size_t nPts = physfield[0].size();
152 Array<OneD, NekDouble> tmp(nPts);
153
154 GetDynamicEnergy(physfield, tmp);
155
156 // Calculate rhoe = E - rho*V^2/2
157 Vmath::Vsub(nPts, physfield[m_spacedim + 1], 1, tmp, 1, energy, 1);
158 // Divide by rho
159 Vmath::Vdiv(nPts, energy, 1, physfield[0], 1, energy, 1);
160}
161
162/**
163 * @brief Compute the specific enthalpy \f$ h = e + p/rho \f$.
164 */
166 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
167 Array<OneD, NekDouble> &enthalpy)
168{
169 size_t nPts = physfield[0].size();
170 Array<OneD, NekDouble> energy(nPts, 0.0);
172
173 GetInternalEnergy(physfield, energy);
174 GetPressure(physfield, pressure);
175
176 // Calculate p/rho
177 Vmath::Vdiv(nPts, pressure, 1, physfield[0], 1, enthalpy, 1);
178 // Calculate h = e + p/rho
179 Vmath::Vadd(nPts, energy, 1, enthalpy, 1, enthalpy, 1);
180}
181
182/**
183 * @brief Compute the velocity field \f$ \mathbf{v} \f$ given the momentum
184 * \f$ \rho\mathbf{v} \f$.
185 *
186 * @param physfield Momentum field.
187 * @param velocity Velocity field.
188 */
190 const Array<OneD, Array<OneD, NekDouble>> &physfield,
192{
193 const size_t nPts = physfield[0].size();
194
195 for (size_t i = 0; i < m_spacedim; ++i)
196 {
197 Vmath::Vdiv(nPts, physfield[1 + i], 1, physfield[0], 1, velocity[i], 1);
198 }
199}
200
201/**
202 * @brief Compute the mach number \f$ M = \| \mathbf{v} \|^2 / c \f$.
203 *
204 * @param physfield Input physical field.
205 * @param soundfield The speed of sound corresponding to physfield.
206 * @param mach The resulting mach number \f$ M \f$.
207 */
209 Array<OneD, NekDouble> &soundspeed,
211{
212 const size_t nPts = physfield[0].size();
213
214 Vmath::Vmul(nPts, physfield[1], 1, physfield[1], 1, mach, 1);
215
216 for (size_t i = 1; i < m_spacedim; ++i)
217 {
218 Vmath::Vvtvp(nPts, physfield[1 + i], 1, physfield[1 + i], 1, mach, 1,
219 mach, 1);
220 }
221
222 Vmath::Vdiv(nPts, mach, 1, physfield[0], 1, mach, 1);
223 Vmath::Vdiv(nPts, mach, 1, physfield[0], 1, mach, 1);
224 Vmath::Vsqrt(nPts, mach, 1, mach, 1);
225
226 Vmath::Vdiv(nPts, mach, 1, soundspeed, 1, mach, 1);
227}
228
229/**
230 * @brief Compute the dynamic viscosity using the Sutherland's law
231 * \f$ \mu = \mu_star * (T / T_star)^3/2 * (1 + C) / (T/T_star + C) \f$,
232 * C : 110. /Tref
233 * Tref : the reference temperature, Tref, should always given in Kelvin,
234 * if non-dimensional should be the reference for non-dimensionalizing
235 * muref : the dynamic viscosity or the 1/Re corresponding to Tref
236 * T_star : m_pInf / (m_rhoInf * m_gasConstant),non-dimensional or dimensional
237 *
238 * WARNING, if this routine is modified the same must be done in the
239 * FieldConvert utility ProcessWSS.cpp (this class should be restructured).
240 *
241 * @param temperature Input temperature.
242 * @param mu The resulting dynamic viscosity.
243 */
246{
247 const size_t nPts = temperature.size();
248
249 for (size_t i = 0; i < nPts; ++i)
250 {
251 mu[i] = GetDynamicViscosity(temperature[i]);
252 }
253}
254
255/**
256 * @brief Compute the dynamic viscosity using the Sutherland's law
257 * \f$ \mu = \mu_star * (T / T_star)^3/2 * (T_star + 110) / (T + 110) \f$,
258 */
260 const Array<OneD, const NekDouble> &temperature,
262{
263 const size_t nPts = temperature.size();
264 NekDouble tmp = 0.0;
265 NekDouble ratio;
266
267 for (size_t i = 0; i < nPts; ++i)
268 {
269 ratio = temperature[i] * m_oneOverT_star;
270 tmp = 0.5 * (ratio + 3.0 * m_TRatioSutherland) /
271 (ratio * (ratio + m_TRatioSutherland));
272 DmuDT[i] = mu[i] * tmp * m_oneOverT_star;
273 }
274}
275
277 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
279{
280 const size_t nPts = physfield[0].size();
281
282 // Getting the velocity vector on the 2D normal space
284
285 Vmath::Zero(Vtot.size(), Vtot, 1);
286
287 for (size_t i = 0; i < m_spacedim; ++i)
288 {
289 velocity[i] = Array<OneD, NekDouble>(nPts);
290 }
291
292 GetVelocityVector(physfield, velocity);
293
294 for (size_t i = 0; i < m_spacedim; ++i)
295 {
296 Vmath::Vvtvp(nPts, velocity[i], 1, velocity[i], 1, Vtot, 1, Vtot, 1);
297 }
298
299 Vmath::Vsqrt(nPts, Vtot, 1, Vtot, 1);
300}
301
302/**
303 * @brief Calculate the pressure using the equation of state.
304 *
305 * @param physfield Input momentum.
306 * @param pressure Computed pressure field.
307 */
309 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
311{
312 size_t nPts = physfield[0].size();
313
314 Array<OneD, NekDouble> energy(nPts);
315 GetInternalEnergy(physfield, energy);
316
317 for (size_t i = 0; i < nPts; ++i)
318 {
319 pressure[i] = m_eos->GetPressure(physfield[0][i], energy[i]);
320 }
321}
322
323/**
324 * @brief Compute the temperature using the equation of state.
325 *
326 * @param physfield Input physical field.
327 * @param temperature The resulting temperature \f$ T \f$.
328 */
330 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
331 Array<OneD, NekDouble> &temperature)
332{
333 size_t nPts = physfield[0].size();
334
335 Array<OneD, NekDouble> energy(nPts);
336 GetInternalEnergy(physfield, energy);
337
338 for (size_t i = 0; i < nPts; ++i)
339 {
340 temperature[i] = m_eos->GetTemperature(physfield[0][i], energy[i]);
341 }
342}
343
344/**
345 * @brief Compute the sound speed using the equation of state.
346 *
347 * @param physfield Input physical field
348 * @param soundspeed The resulting sound speed \f$ c \f$.
349 */
351 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
352 Array<OneD, NekDouble> &soundspeed)
353{
354 size_t nPts = physfield[0].size();
355
356 Array<OneD, NekDouble> energy(nPts);
357 GetInternalEnergy(physfield, energy);
358
359 for (size_t i = 0; i < nPts; ++i)
360 {
361 soundspeed[i] = m_eos->GetSoundSpeed(physfield[0][i], energy[i]);
362 }
363}
364
365/**
366 * @brief Compute the entropy using the equation of state.
367 *
368 * @param physfield Input physical field
369 * @param soundspeed The resulting sound speed \f$ c \f$.
370 */
372 const Array<OneD, const Array<OneD, NekDouble>> &physfield,
373 Array<OneD, NekDouble> &entropy)
374{
375 size_t nPts = physfield[0].size();
376
377 Array<OneD, NekDouble> energy(nPts);
378 GetInternalEnergy(physfield, energy);
379
380 for (size_t i = 0; i < nPts; ++i)
381 {
382 entropy[i] = m_eos->GetEntropy(physfield[0][i], energy[i]);
383 }
384}
385
386/**
387 * @brief Compute \f$ e(rho,p) \f$ using the equation of state.
388 *
389 * @param rho Input density
390 * @param pressure Input pressure
391 * @param energy The resulting internal energy.
392 */
396{
397 size_t nPts = rho.size();
398
399 for (size_t i = 0; i < nPts; ++i)
400 {
401 energy[i] = m_eos->GetEFromRhoP(rho[i], pressure[i]);
402 }
403}
404
405/**
406 * @brief Compute \f$ rho(p,T) \f$ using the equation of state.
407 *
408 * @param pressure Input pressure
409 * @param temperature Input temperature
410 * @param rho The resulting density
411 */
413 const Array<OneD, NekDouble> &temperature,
415{
416 size_t nPts = pressure.size();
417
418 for (size_t i = 0; i < nPts; ++i)
419 {
420 rho[i] = m_eos->GetRhoFromPT(pressure[i], temperature[i]);
421 }
422}
423
426 const Array<OneD, const Array<OneD, NekDouble>> &consVar,
427 const Array<OneD, NekDouble> &div,
428 const Array<OneD, NekDouble> &curlSquared)
429{
430 size_t nTracePts = fields[0]->GetTrace()->GetTotPoints();
432 {
433 size_t nPts = fields[0]->GetTotPoints();
434 m_muAv = Array<OneD, NekDouble>(nPts, 0.0);
435 m_muAvTrace = Array<OneD, NekDouble>(nTracePts, 0.0);
436 SetElmtMinHP(fields);
437 }
438
439 if (m_shockSensorType == "Modal")
440 {
441 // Get viscosity based on modal sensor
442 GetMuAv(fields, consVar, m_muAv);
443 }
444 else
445 {
446 // Get viscosity based on dilatation sensor
447 GetMuAv(fields, consVar, div, m_muAv);
448 }
449
450 // Apply Ducros sensor
451 if (m_ducrosSensor != "Off")
452 {
453 ApplyDucros(div, curlSquared, m_muAv);
454 }
455
456 // Apply approximate C0 smoothing
457 if (m_smoothing == "C0")
458 {
460 }
461
462 // Set trace AV
463 Array<OneD, NekDouble> muFwd(nTracePts, 0.0), muBwd(nTracePts, 0.0);
464 fields[0]->GetFwdBwdTracePhys(m_muAv, muFwd, muBwd, false, false, false);
465 for (size_t p = 0; p < nTracePts; ++p)
466 {
467 m_muAvTrace[p] = 0.5 * (muFwd[p] + muBwd[p]);
468 }
469}
470
471/**
472 * @brief Compute an estimate of minimum h/p
473 * for each element of the expansion.
474 */
477{
478 size_t nElements = fields[0]->GetExpSize();
480 {
481 m_hOverP = Array<OneD, NekDouble>(nElements, 1.0);
482 }
483
484 // Determine h/p scaling
485 Array<OneD, int> pOrderElmt = fields[0]->EvalBasisNumModesMaxPerExp();
486 int expdim = fields[0]->GetGraph()->GetMeshDimension();
487 for (size_t e = 0; e < nElements; e++)
488 {
489 NekDouble h = 1.0e+10;
490 switch (expdim)
491 {
492 case 3:
493 {
495 exp3D = fields[0]->GetExp(e)->as<LocalRegions::Expansion3D>();
496 for (int i = 0; i < exp3D->GetNtraces(); ++i)
497 {
498 h = std::min(
499 h, exp3D->GetGeom3D()->GetEdge(i)->GetVertex(0)->dist(*(
500 exp3D->GetGeom3D()->GetEdge(i)->GetVertex(1))));
501 }
502 break;
503 }
504
505 case 2:
506 {
508 exp2D = fields[0]->GetExp(e)->as<LocalRegions::Expansion2D>();
509 for (int i = 0; i < exp2D->GetNtraces(); ++i)
510 {
511 h = std::min(
512 h, exp2D->GetGeom2D()->GetEdge(i)->GetVertex(0)->dist(*(
513 exp2D->GetGeom2D()->GetEdge(i)->GetVertex(1))));
514 }
515 break;
516 }
517 case 1:
518 {
520 exp1D = fields[0]->GetExp(e)->as<LocalRegions::Expansion1D>();
521
522 h = std::min(h, exp1D->GetGeom1D()->GetVertex(0)->dist(
523 *(exp1D->GetGeom1D()->GetVertex(1))));
524
525 break;
526 }
527 default:
528 {
529 ASSERTL0(false, "Dimension out of bound.")
530 }
531 }
532
533 // Store h/p scaling
534 m_hOverP[e] = h / std::max(pOrderElmt[e] - 1, 1);
535 }
536}
537
539{
540 ASSERTL1(m_hOverP != NullNekDouble1DArray, "m_hOverP not set");
541 return m_hOverP;
542}
543
546 const Array<OneD, const Array<OneD, NekDouble>> &physarray,
548 int offset)
549{
550 NekDouble Skappa;
551 NekDouble order;
553 Array<OneD, int> expOrderElement = field->EvalBasisNumModesMaxPerExp();
554
555 for (int e = 0; e < field->GetExpSize(); e++)
556 {
557 int numModesElement = expOrderElement[e];
558 int nElmtPoints = field->GetExp(e)->GetTotPoints();
559 int physOffset = field->GetPhys_Offset(e);
560 int nElmtCoeffs = field->GetExp(e)->GetNcoeffs();
561 int numCutOff = numModesElement - offset;
562
563 if (numModesElement <= offset)
564 {
565 Vmath::Fill(nElmtPoints, 0.0, tmp = Sensor + physOffset, 1);
566 Vmath::Fill(nElmtPoints, 0.0, tmp = SensorKappa + physOffset, 1);
567 continue;
568 }
569
570 // create vector to save the solution points per element at P = p;
571 Array<OneD, NekDouble> elmtPhys(nElmtPoints,
572 tmp = physarray[0] + physOffset);
573 // Compute coefficients
574 Array<OneD, NekDouble> elmtCoeffs(nElmtCoeffs, 0.0);
575 field->GetExp(e)->FwdTrans(elmtPhys, elmtCoeffs);
576
577 // ReduceOrderCoeffs reduces the polynomial order of the solution
578 // that is represented by the coeffs given as an inarray. This is
579 // done by projecting the higher order solution onto the orthogonal
580 // basis and padding the higher order coefficients with zeros.
581 Array<OneD, NekDouble> reducedElmtCoeffs(nElmtCoeffs, 0.0);
582 field->GetExp(e)->ReduceOrderCoeffs(numCutOff, elmtCoeffs,
583 reducedElmtCoeffs);
584
585 Array<OneD, NekDouble> reducedElmtPhys(nElmtPoints, 0.0);
586 field->GetExp(e)->BwdTrans(reducedElmtCoeffs, reducedElmtPhys);
587
588 NekDouble numerator = 0.0;
589 NekDouble denominator = 0.0;
590
591 // Determining the norm of the numerator of the Sensor
592 Array<OneD, NekDouble> difference(nElmtPoints, 0.0);
593 Vmath::Vsub(nElmtPoints, elmtPhys, 1, reducedElmtPhys, 1, difference,
594 1);
595
596 numerator = Vmath::Dot(nElmtPoints, difference, difference);
597 denominator = Vmath::Dot(nElmtPoints, elmtPhys, elmtPhys);
598
599 NekDouble elmtSensor = sqrt(numerator / denominator);
600 elmtSensor =
601 log10(std::max(elmtSensor, NekConstants::kNekMachineEpsilon));
602
603 Vmath::Fill(nElmtPoints, elmtSensor, tmp = Sensor + physOffset, 1);
604
605 // Compute reference value for sensor
606 order = std::max(numModesElement - 1, 1);
607 if (order > 0)
608 {
609 Skappa = m_Skappa - 4.25 * log10(static_cast<NekDouble>(order));
610 }
611 else
612 {
613 Skappa = 0.0;
614 }
615
616 // Compute artificial viscosity
617 NekDouble elmtSensorKappa;
618 if (elmtSensor < (Skappa - m_Kappa))
619 {
620 elmtSensorKappa = 0;
621 }
622 else if (elmtSensor > (Skappa + m_Kappa))
623 {
624 elmtSensorKappa = 1.0;
625 }
626 else
627 {
628 elmtSensorKappa =
629 0.5 * (1 + sin(M_PI * (elmtSensor - Skappa) / (2 * m_Kappa)));
630 }
631 Vmath::Fill(nElmtPoints, elmtSensorKappa,
632 tmp = SensorKappa + physOffset, 1);
633 }
634}
635
636/**
637 * @brief Calculate the physical artificial viscosity based on modal sensor.
638 *
639 * @param consVar Input field.
640 */
643 const Array<OneD, const Array<OneD, NekDouble>> &consVar,
645{
646 size_t nPts = consVar[0].size();
647 // Determine the maximum wavespeed
648 Array<OneD, NekDouble> Lambdas(nPts, 0.0);
649 Array<OneD, NekDouble> soundspeed(nPts, 0.0);
650 Array<OneD, NekDouble> absVelocity(nPts, 0.0);
651 GetSoundSpeed(consVar, soundspeed);
652 GetAbsoluteVelocity(consVar, absVelocity);
653 Vmath::Vadd(nPts, absVelocity, 1, soundspeed, 1, Lambdas, 1);
654
655 // Compute sensor based on rho
656 Array<OneD, NekDouble> Sensor(nPts, 0.0);
657 GetSensor(fields[0], consVar, Sensor, muAv, 1);
658
660 size_t nElmt = fields[0]->GetExpSize();
661 for (size_t e = 0; e < nElmt; ++e)
662 {
663 int physOffset = fields[0]->GetPhys_Offset(e);
664 int nElmtPoints = fields[0]->GetExp(e)->GetTotPoints();
665
666 // Compute the maximum wave speed
667 NekDouble LambdaElmt = 0.0;
668 LambdaElmt = Vmath::Vmax(nElmtPoints, tmp = Lambdas + physOffset, 1);
669
670 // Compute average bounded density
671 NekDouble rhoAve =
672 Vmath::Vsum(nElmtPoints, tmp = consVar[0] + physOffset, 1);
673 rhoAve = rhoAve / nElmtPoints;
674 rhoAve = Smath::Smax(rhoAve, 1.0e-4, 1.0e+4);
675
676 // Scale sensor by coeff, h/p, and density
677 LambdaElmt *= m_mu0 * m_hOverP[e] * rhoAve;
678 Vmath::Smul(nElmtPoints, LambdaElmt, tmp = muAv + physOffset, 1,
679 tmp = muAv + physOffset, 1);
680 }
681}
682
683/**
684 * @brief Calculate the physical artificial viscosity based on dilatation of
685 * velocity vector.
686 *
687 * @param
688 */
691 const Array<OneD, const Array<OneD, NekDouble>> &consVar,
693{
694 size_t nPts = consVar[0].size();
695
696 // Get sound speed
697 // theoretically it should be used the critical sound speed, this
698 // matters for large Mach numbers (above 3.0)
699 Array<OneD, NekDouble> soundSpeed(nPts, 0.0);
700 GetSoundSpeed(consVar, soundSpeed);
701
702 // Get abosolute velocity to compute lambda
703 Array<OneD, NekDouble> absVelocity(nPts, 0.0);
704 GetAbsoluteVelocity(consVar, absVelocity);
705
706 // Loop over elements
707 size_t nElmt = fields[0]->GetExpSize();
708 for (size_t e = 0; e < nElmt; ++e)
709 {
710 int nElmtPoints = fields[0]->GetExp(e)->GetTotPoints();
711 int physOffset = fields[0]->GetPhys_Offset(e);
712 int physEnd = physOffset + nElmtPoints;
713
714 NekDouble hOpTmp = m_hOverP[e];
715
716 // Loop over the points
717 for (int p = physOffset; p < physEnd; ++p)
718 {
719 // Get non-dimensional sensor based on dilatation
720 NekDouble sSpeedTmp = soundSpeed[p];
721 // (only compression waves)
722 NekDouble divTmp = -div[p];
723 divTmp = std::max(divTmp, 0.0);
724 NekDouble sensor = m_mu0 * hOpTmp * divTmp / sSpeedTmp;
725 // Scale to viscosity scale
726 NekDouble rho = consVar[0][p];
727 NekDouble lambda = sSpeedTmp + absVelocity[p];
728 muAv[p] = sensor * rho * lambda * hOpTmp;
729 }
730 }
731}
732
733/**
734 * @brief Apply Ducros (anti-vorticity) sensor averaged over the element.
735 *
736 * @param field Input Field
737 */
739 const Array<OneD, NekDouble> &curlSquare,
741{
742 // machine eps**2
743 NekDouble eps = std::numeric_limits<NekDouble>::epsilon();
744 eps *= eps;
745
746 // loop over points
747 size_t nPts = div.size();
748 for (size_t p = 0; p < nPts; ++p)
749 {
750 NekDouble tmpDiv2 = div[p];
751 tmpDiv2 *= tmpDiv2;
752 NekDouble denDuc = tmpDiv2 + curlSquare[p] + eps;
753 NekDouble Duc = tmpDiv2 / denDuc;
754 // apply
755 muAv[p] *= Duc;
756 }
757}
758
759/**
760 * @brief Make field C0.
761 *
762 * @param field Input Field
763 */
765{
766 // Make sure that the C0 projection operator has been allocated. Note that
767 // the VariableConverter object can be allocated without the C0 smoother.
768 // This is why this check is needed. Ideally, the C0 smoother is separated
769 // from the VariableConverter class.
770 ASSERTL0(m_C0ProjectExp.get() != nullptr,
771 "C0 projection operator not initialized in "
772 "VariableConverter::ApplyC0Smooth()");
773
774 int nCoeffs = m_C0ProjectExp->GetNcoeffs();
775 Array<OneD, NekDouble> muFwd(nCoeffs);
776 Array<OneD, NekDouble> weights(nCoeffs, 1.0);
777 // Assemble global expansion coefficients for viscosity
778 m_C0ProjectExp->FwdTransLocalElmt(field, m_C0ProjectExp->UpdateCoeffs());
779 m_C0ProjectExp->Assemble();
780 Vmath::Vcopy(nCoeffs, m_C0ProjectExp->GetCoeffs(), 1, muFwd, 1);
781 // Global coefficients
782 Vmath::Vcopy(nCoeffs, weights, 1, m_C0ProjectExp->UpdateCoeffs(), 1);
783 // This is the sign vector
784 m_C0ProjectExp->GlobalToLocal();
785 // Get weights
786 m_C0ProjectExp->Assemble();
787 // Divide
788 Vmath::Vdiv(nCoeffs, muFwd, 1, m_C0ProjectExp->GetCoeffs(), 1,
789 m_C0ProjectExp->UpdateCoeffs(), 1);
790 // Get local coefficients
791 m_C0ProjectExp->GlobalToLocal();
792 // Get C0 field
793 m_C0ProjectExp->BwdTrans(m_C0ProjectExp->GetCoeffs(), field);
794}
795
796} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define WARNINGL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
SpatialDomains::Geometry1DSharedPtr GetGeom1D() const
Definition: Expansion1D.h:104
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void ApplyC0Smooth(Array< OneD, NekDouble > &field)
Make field C0.
void GetInternalEnergy(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &energy)
Compute the specific internal energy .
EquationOfStateSharedPtr m_eos
void GetSensor(const MultiRegions::ExpListSharedPtr &field, const Array< OneD, const Array< OneD, NekDouble > > &physarray, Array< OneD, NekDouble > &Sensor, Array< OneD, NekDouble > &SensorKappa, int offset=1)
Array< OneD, NekDouble > m_hOverP
h/p scaling
void GetRhoFromPT(const Array< OneD, NekDouble > &pressure, const Array< OneD, NekDouble > &temperature, Array< OneD, NekDouble > &rho)
Compute using the equation of state.
void GetAbsoluteVelocity(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &Vtot)
void GetMuAv(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &muAv)
Calculate the physical artificial viscosity based on modal sensor.
void SetAv(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const Array< OneD, NekDouble > > &consVar, const Array< OneD, NekDouble > &div=NullNekDouble1DArray, const Array< OneD, NekDouble > &curlSquared=NullNekDouble1DArray)
Array< OneD, NekDouble > m_muAvTrace
void GetDmuDT(const Array< OneD, const NekDouble > &temperature, const Array< OneD, const NekDouble > &mu, Array< OneD, NekDouble > &DmuDT)
Compute the dynamic viscosity using the Sutherland's law ,.
void GetEFromRhoP(const Array< OneD, NekDouble > &rho, const Array< OneD, NekDouble > &pressure, Array< OneD, NekDouble > &energy)
Compute using the equation of state.
void SetElmtMinHP(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields)
Compute an estimate of minimum h/p for each element of the expansion.
MultiRegions::ContFieldSharedPtr m_C0ProjectExp
void GetPressure(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &pressure)
Calculate the pressure using the equation of state.
LibUtilities::SessionReaderSharedPtr m_session
NekDouble m_mu0
Shock sensor.
void GetTemperature(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &temperature)
Compute the temperature using the equation of state.
void ApplyDucros(const Array< OneD, NekDouble > &div, const Array< OneD, NekDouble > &curlSquare, Array< OneD, NekDouble > &muAv)
Apply Ducros (anti-vorticity) sensor averaged over the element.
void GetDynamicEnergy(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &energy)
Compute the dynamic energy .
void GetSoundSpeed(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &soundspeed)
Compute the sound speed using the equation of state.
VariableConverter(const LibUtilities::SessionReaderSharedPtr &pSession, const int spaceDim, const SpatialDomains::MeshGraphSharedPtr &pGraph=nullptr)
void GetEntropy(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &entropy)
Compute the entropy using the equation of state.
void GetMach(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &soundspeed, Array< OneD, NekDouble > &mach)
Compute the mach number .
void GetVelocityVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &velocity)
Compute the velocity field given the momentum .
void GetEnthalpy(const Array< OneD, const Array< OneD, NekDouble > > &physfield, Array< OneD, NekDouble > &enthalpy)
Compute the specific enthalpy .
Array< OneD, NekDouble > & GetElmtMinHP()
Array< OneD, NekDouble > m_muAv
storage
void GetDynamicViscosity(const Array< OneD, const NekDouble > &temperature, Array< OneD, NekDouble > &mu)
Compute the dynamic viscosity using the Sutherland's law , C : 110. /Tref Tref : the reference temper...
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:46
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:50
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:47
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
static const NekDouble kNekMachineEpsilon
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
EquationOfStateFactory & GetEquationOfStateFactory()
Declaration of the equation of state factory singleton.
static Array< OneD, NekDouble > NullNekDouble1DArray
double NekDouble
T Smax(const T a, const T b, const T k)
Return the soft max of between two scalars.
Definition: Smath.hpp:53
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
T Dot(int n, const T *w, const T *x)
dot product
Definition: Vmath.hpp:761
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.hpp:180
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.hpp:126
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.hpp:54
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.hpp:644
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.hpp:220
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285