Nektar++
VCSMapping.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: VCSMapping.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: Velocity Correction Scheme w/ coordinate transformation
32// for the Incompressible Navier Stokes equations
33///////////////////////////////////////////////////////////////////////////////
34
38
39#include <boost/algorithm/string.hpp>
40
41namespace Nektar
42{
43
44std::string VCSMapping::className =
46 "VCSMapping", VCSMapping::create);
47
49 LibUtilities::SessionReader::RegisterEnumValue("SolverType", "VCSMapping",
51
52/**
53 * Constructor. Creates ...
54 *
55 * \param
56 * \param
57 */
60 : UnsteadySystem(pSession, pGraph),
61 VelocityCorrectionScheme(pSession, pGraph)
62{
63}
64
65void VCSMapping::v_InitObject(bool DeclareField)
66{
68
70 ASSERTL0(m_mapping, "Could not create mapping in VCSMapping.");
71
72 std::string vExtrapolation = "Mapping";
74 vExtrapolation, m_session, m_fields, m_pressure, m_velocity,
76 m_extrapolation->SubSteppingTimeIntegration(m_intScheme);
77 m_extrapolation->GenerateBndElmtExpansion();
78 m_extrapolation->GenerateHOPBCMap(m_session);
79
80 // Storage to extrapolate pressure forcing
81 size_t physTot = m_fields[0]->GetTotPoints();
82 size_t intSteps = 1;
83
84 if (m_intScheme->GetName() == "IMEX" ||
85 m_intScheme->GetName() == "IMEXGear")
86 {
87 m_intSteps = m_intScheme->GetOrder();
88 }
89 else
90 {
92 "Integration method not suitable: "
93 "Options include IMEXGear or IMEXOrder{1,2,3,4}");
94 }
95
97 for (size_t i = 0; i < m_presForcingCorrection.size(); i++)
98 {
100 }
101 m_verbose = (m_session->DefinesCmdLineArgument("verbose")) ? true : false;
102
103 // Load solve parameters related to the mapping
104 // Flags determining if pressure/viscous terms should be treated implicitly
105 m_session->MatchSolverInfo("MappingImplicitPressure", "True",
106 m_implicitPressure, false);
107 m_session->MatchSolverInfo("MappingImplicitViscous", "True",
108 m_implicitViscous, false);
109 m_session->MatchSolverInfo("MappingNeglectViscous", "True",
110 m_neglectViscous, false);
111
113 {
114 m_implicitViscous = false;
115 }
116
117 // Tolerances and relaxation parameters for implicit terms
118 m_session->LoadParameter("MappingPressureTolerance", m_pressureTolerance,
119 1e-12);
120 m_session->LoadParameter("MappingViscousTolerance", m_viscousTolerance,
121 1e-12);
122 m_session->LoadParameter("MappingPressureRelaxation", m_pressureRelaxation,
123 1.0);
124 m_session->LoadParameter("MappingViscousRelaxation", m_viscousRelaxation,
125 1.0);
126}
127
128void VCSMapping::v_DoInitialise(bool dumpInitialConditions)
129{
130 UnsteadySystem::v_DoInitialise(dumpInitialConditions);
131
132 // Set up Field Meta Data for output files
133 m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
134 m_fieldMetaDataMap["TimeStep"] =
135 boost::lexical_cast<std::string>(m_timestep);
136
137 // Correct Dirichlet boundary conditions to account for mapping
138 m_mapping->UpdateBCs(0.0);
139 //
141 for (int i = 0; i < m_nConvectiveFields; ++i)
142 {
143 m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
144 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
145 m_fields[i]->UpdatePhys());
147 }
148
149 // Initialise m_gradP
150 size_t physTot = m_fields[0]->GetTotPoints();
152 for (int i = 0; i < m_nConvectiveFields; ++i)
153 {
154 m_gradP[i] = Array<OneD, NekDouble>(physTot, 0.0);
156 m_pressure->GetPhys(), m_gradP[i]);
157 if (m_pressure->GetWaveSpace())
158 {
159 m_pressure->HomogeneousBwdTrans(physTot, m_gradP[i], m_gradP[i]);
160 }
161 }
162}
163
164/**
165 * Explicit part of the method - Advection, Forcing + HOPBCs
166 */
168 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
169 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
170{
171 EvaluateAdvectionTerms(inarray, outarray, time);
172
173 // Smooth advection
175 {
176 for (int i = 0; i < m_nConvectiveFields; ++i)
177 {
178 m_pressure->SmoothField(outarray[i]);
179 }
180 }
181
182 // Add forcing terms
183 for (auto &x : m_forcing)
184 {
185 x->Apply(m_fields, inarray, outarray, time);
186 }
187
188 // Add mapping terms
189 ApplyIncNSMappingForcing(inarray, outarray);
190
191 // Calculate High-Order pressure boundary conditions
192 m_extrapolation->EvaluatePressureBCs(inarray, outarray, m_kinvis);
193
194 // Update mapping and deal with Dirichlet boundary conditions
195 if (m_mapping->IsTimeDependent())
196 {
197 if (m_mapping->IsFromFunction())
198 {
199 // If the transformation is explicitly defined, update it here
200 // Otherwise, it will be done somewhere else (ForcingMovingBody)
201 m_mapping->UpdateMapping(time + m_timestep);
202 }
203 m_mapping->UpdateBCs(time + m_timestep);
204 }
205}
206
207/**
208 * Forcing term for Poisson solver solver
209 */
211 const Array<OneD, const Array<OneD, NekDouble>> &fields,
213{
214 if (m_mapping->HasConstantJacobian())
215 {
217 aii_Dt);
218 }
219 else
220 {
221 size_t physTot = m_fields[0]->GetTotPoints();
222 size_t nvel = m_nConvectiveFields;
223 Array<OneD, NekDouble> wk(physTot, 0.0);
224
225 Array<OneD, NekDouble> Jac(physTot, 0.0);
226 m_mapping->GetJacobian(Jac);
227
228 // Calculate div(J*u/Dt)
229 Vmath::Zero(physTot, Forcing[0], 1);
230 for (size_t i = 0; i < nvel; ++i)
231 {
232 if (m_fields[i]->GetWaveSpace())
233 {
234 m_fields[i]->HomogeneousBwdTrans(physTot, fields[i], wk);
235 }
236 else
237 {
238 Vmath::Vcopy(physTot, fields[i], 1, wk, 1);
239 }
240 Vmath::Vmul(physTot, wk, 1, Jac, 1, wk, 1);
241 if (m_fields[i]->GetWaveSpace())
242 {
243 m_fields[i]->HomogeneousFwdTrans(physTot, wk, wk);
244 }
245 m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[i], wk, wk);
246 Vmath::Vadd(physTot, wk, 1, Forcing[0], 1, Forcing[0], 1);
247 }
248 Vmath::Smul(physTot, 1.0 / aii_Dt, Forcing[0], 1, Forcing[0], 1);
249
250 //
251 // If the mapping viscous terms are being treated explicitly
252 // we need to apply a correction to the forcing
254 {
255 bool wavespace = m_fields[0]->GetWaveSpace();
256 m_fields[0]->SetWaveSpace(false);
257
258 //
259 // Part 1: div(J*grad(U/J . grad(J)))
262 for (size_t i = 0; i < tmp.size(); i++)
263 {
264 tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
265 velocity[i] = Array<OneD, NekDouble>(physTot, 0.0);
266 if (wavespace)
267 {
268 m_fields[0]->HomogeneousBwdTrans(
269 physTot, m_fields[i]->GetPhys(), velocity[i]);
270 }
271 else
272 {
273 Vmath::Vcopy(physTot, m_fields[i]->GetPhys(), 1,
274 velocity[i], 1);
275 }
276 }
277 // Calculate wk = U.grad(J)
278 m_mapping->DotGradJacobian(velocity, wk);
279 // Calculate wk = (U.grad(J))/J
280 Vmath::Vdiv(physTot, wk, 1, Jac, 1, wk, 1);
281 // J*grad[(U.grad(J))/J]
282 for (size_t i = 0; i < nvel; ++i)
283 {
284 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], wk,
285 tmp[i]);
286 Vmath::Vmul(physTot, Jac, 1, tmp[i], 1, tmp[i], 1);
287 }
288 // div(J*grad[(U.grad(J))/J])
289 Vmath::Zero(physTot, wk, 1);
290 for (size_t i = 0; i < nvel; ++i)
291 {
292 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], tmp[i],
293 tmp[i]);
294 Vmath::Vadd(physTot, wk, 1, tmp[i], 1, wk, 1);
295 }
296
297 // Part 2: grad(J) . curl(curl(U))
298 m_mapping->CurlCurlField(velocity, tmp, m_implicitViscous);
299 // dont need velocity any more, so reuse it
300 m_mapping->DotGradJacobian(tmp, velocity[0]);
301
302 // Add two parts
303 Vmath::Vadd(physTot, velocity[0], 1, wk, 1, wk, 1);
304
305 // Multiply by kinvis and prepare to extrapolate
306 size_t nlevels = m_presForcingCorrection.size();
307 Vmath::Smul(physTot, m_kinvis, wk, 1,
308 m_presForcingCorrection[nlevels - 1], 1);
309
310 // Extrapolate correction
311 m_extrapolation->ExtrapolateArray(m_presForcingCorrection);
312
313 // Put in wavespace
314 if (wavespace)
315 {
316 m_fields[0]->HomogeneousFwdTrans(
317 physTot, m_presForcingCorrection[nlevels - 1], wk);
318 }
319 else
320 {
321 Vmath::Vcopy(physTot, m_presForcingCorrection[nlevels - 1], 1,
322 wk, 1);
323 }
324 // Apply correction: Forcing = Forcing - correction
325 Vmath::Vsub(physTot, Forcing[0], 1, wk, 1, Forcing[0], 1);
326
327 m_fields[0]->SetWaveSpace(wavespace);
328 }
329 }
330}
331
332/**
333 * Forcing term for Helmholtz solver
334 */
336 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
338{
339 NekDouble aii_dtinv = 1.0 / aii_Dt;
340 size_t physTot = m_fields[0]->GetTotPoints();
341
342 // Grad p
343 m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
344
345 size_t nvel = m_velocity.size();
346 if (nvel == 2)
347 {
348 m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1]);
349 }
350 else
351 {
352 m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1],
353 Forcing[2]);
354 }
355
356 // Copy grad p in physical space to m_gradP to reuse later
357 if (m_pressure->GetWaveSpace())
358 {
359 for (size_t i = 0; i < nvel; i++)
360 {
361 m_pressure->HomogeneousBwdTrans(physTot, Forcing[i], m_gradP[i]);
362 }
363 }
364 else
365 {
366 for (size_t i = 0; i < nvel; i++)
367 {
368 Vmath::Vcopy(physTot, Forcing[i], 1, m_gradP[i], 1);
369 }
370 }
371
372 if ((!m_mapping->HasConstantJacobian()) || m_implicitPressure)
373 {
374 // If pressure terms are treated explicitly, we need to divide by J
375 // if they are implicit, we need to calculate G(p)
377 {
378 m_mapping->RaiseIndex(m_gradP, Forcing);
379 }
380 else
381 {
382 Array<OneD, NekDouble> Jac(physTot, 0.0);
383 m_mapping->GetJacobian(Jac);
384 for (size_t i = 0; i < nvel; i++)
385 {
386 Vmath::Vdiv(physTot, m_gradP[i], 1, Jac, 1, Forcing[i], 1);
387 }
388 }
389 // Transform back to wavespace
390 if (m_pressure->GetWaveSpace())
391 {
392 for (size_t i = 0; i < nvel; i++)
393 {
394 m_pressure->HomogeneousFwdTrans(physTot, Forcing[i],
395 Forcing[i]);
396 }
397 }
398 }
399
400 // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
401 // need to be updated for the convected fields.
402 for (size_t i = 0; i < nvel; ++i)
403 {
404 Blas::Daxpy(physTot, -aii_dtinv, inarray[i], 1, Forcing[i], 1);
405 Blas::Dscal(physTot, 1.0 / m_kinvis, &(Forcing[i])[0], 1);
406 }
407}
408
409/**
410 * Solve pressure system
411 */
413{
415 {
417 }
418 else
419 {
420 size_t physTot = m_fields[0]->GetTotPoints();
421 size_t nvel = m_nConvectiveFields;
422 bool converged = false; // flag to mark if system converged
423 size_t s = 0; // iteration counter
424 NekDouble error; // L2 error at current iteration
425 NekDouble forcing_L2 = 0.0; // L2 norm of F
426
427 size_t maxIter;
428 m_session->LoadParameter("MappingMaxIter", maxIter, 5000);
429
430 // rhs of the equation at current iteration
431 Array<OneD, NekDouble> F_corrected(physTot, 0.0);
432 // Pressure field at previous iteration
433 Array<OneD, NekDouble> previous_iter(physTot, 0.0);
434 // Temporary variables
438 for (size_t i = 0; i < nvel; ++i)
439 {
440 wk1[i] = Array<OneD, NekDouble>(physTot, 0.0);
441 wk2[i] = Array<OneD, NekDouble>(physTot, 0.0);
442 gradP[i] = Array<OneD, NekDouble>(physTot, 0.0);
443 }
444
445 // Jacobian
446 Array<OneD, NekDouble> Jac(physTot, 0.0);
447 m_mapping->GetJacobian(Jac);
448
449 // Factors for Laplacian system
452
453 m_pressure->BwdTrans(m_pressure->GetCoeffs(), m_pressure->UpdatePhys());
454 forcing_L2 = m_pressure->L2(Forcing, wk1[0]);
455 while (!converged)
456 {
457 // Update iteration counter and set previous iteration field
458 // (use previous timestep solution for first iteration)
459 s++;
460 ASSERTL0(s < maxIter,
461 "VCSMapping exceeded maximum number of iterations.");
462
463 Vmath::Vcopy(physTot, m_pressure->GetPhys(), 1, previous_iter, 1);
464
465 // Correct pressure bc to account for iteration
466 m_extrapolation->CorrectPressureBCs(previous_iter);
467
468 //
469 // Calculate forcing term for this iteration
470 //
471 for (size_t i = 0; i < nvel; ++i)
472 {
474 previous_iter, gradP[i]);
475 if (m_pressure->GetWaveSpace())
476 {
477 m_pressure->HomogeneousBwdTrans(physTot, gradP[i], wk1[i]);
478 }
479 else
480 {
481 Vmath::Vcopy(physTot, gradP[i], 1, wk1[i], 1);
482 }
483 }
484 m_mapping->RaiseIndex(wk1, wk2); // G(p)
485
486 m_mapping->Divergence(wk2, F_corrected); // div(G(p))
487 if (!m_mapping->HasConstantJacobian())
488 {
489 Vmath::Vmul(physTot, F_corrected, 1, Jac, 1, F_corrected, 1);
490 }
491 // alpha*J*div(G(p))
492 Vmath::Smul(physTot, m_pressureRelaxation, F_corrected, 1,
493 F_corrected, 1);
494 if (m_pressure->GetWaveSpace())
495 {
496 m_pressure->HomogeneousFwdTrans(physTot, F_corrected,
497 F_corrected);
498 }
499 // alpha*J*div(G(p)) - p_ii
500 for (int i = 0; i < m_nConvectiveFields; ++i)
501 {
503 gradP[i], wk1[0]);
504 Vmath::Vsub(physTot, F_corrected, 1, wk1[0], 1, F_corrected, 1);
505 }
506 // p_i,i - J*div(G(p))
507 Vmath::Neg(physTot, F_corrected, 1);
508 // alpha*F - alpha*J*div(G(p)) + p_i,i
509 Vmath::Smul(physTot, m_pressureRelaxation, Forcing, 1, wk1[0], 1);
510 Vmath::Vadd(physTot, wk1[0], 1, F_corrected, 1, F_corrected, 1);
511
512 //
513 // Solve system
514 //
515 m_pressure->HelmSolve(F_corrected, m_pressure->UpdateCoeffs(),
516 factors);
517 m_pressure->BwdTrans(m_pressure->GetCoeffs(),
518 m_pressure->UpdatePhys());
519
520 //
521 // Test convergence
522 //
523 error = m_pressure->L2(m_pressure->GetPhys(), previous_iter);
524 if (forcing_L2 != 0)
525 {
526 if ((error / forcing_L2 < m_pressureTolerance))
527 {
528 converged = true;
529 }
530 }
531 else
532 {
533 if (error < m_pressureTolerance)
534 {
535 converged = true;
536 }
537 }
538 }
539 if (m_verbose && m_session->GetComm()->GetRank() == 0)
540 {
541 std::cout << " Pressure system (mapping) converged in " << s
542 << " iterations with error = " << error << std::endl;
543 }
544 }
545}
546
547/**
548 * Solve velocity system
549 */
552 [[maybe_unused]] const Array<OneD, const Array<OneD, NekDouble>> &inarray,
553 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble aii_Dt)
554{
556 {
558 aii_Dt);
559 }
560 else
561 {
562 size_t physTot = m_fields[0]->GetTotPoints();
563 size_t nvel = m_nConvectiveFields;
564 bool converged = false; // flag to mark if system converged
565 size_t s = 0; // iteration counter
566 NekDouble error, max_error; // L2 error at current iteration
567
568 size_t maxIter;
569 m_session->LoadParameter("MappingMaxIter", maxIter, 5000);
570
571 // L2 norm of F
573
574 // rhs of the equation at current iteration
575 Array<OneD, Array<OneD, NekDouble>> F_corrected(nvel);
576 // Solution at previous iteration
577 Array<OneD, Array<OneD, NekDouble>> previous_iter(nvel);
578 // Working space
580 for (size_t i = 0; i < nvel; ++i)
581 {
582 F_corrected[i] = Array<OneD, NekDouble>(physTot, 0.0);
583 previous_iter[i] = Array<OneD, NekDouble>(physTot, 0.0);
584 wk[i] = Array<OneD, NekDouble>(physTot, 0.0);
585 }
586
587 // Factors for Helmholtz system
590 1.0 * m_viscousRelaxation / aii_Dt / m_kinvis;
592 {
596 }
597
598 // Calculate L2-norm of F and set initial solution for iteration
599 for (size_t i = 0; i < nvel; ++i)
600 {
601 forcing_L2[i] = m_fields[0]->L2(Forcing[i], wk[0]);
602 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), previous_iter[i]);
603 }
604
605 while (!converged)
606 {
607 converged = true;
608 // Iteration counter
609 s++;
610 ASSERTL0(s < maxIter,
611 "VCSMapping exceeded maximum number of iterations.");
612
613 max_error = 0.0;
614
615 //
616 // Calculate forcing term for next iteration
617 //
618
619 // Calculate L(U)- in this parts all components might be coupled
620 if (m_fields[0]->GetWaveSpace())
621 {
622 for (size_t i = 0; i < nvel; ++i)
623 {
624 m_fields[0]->HomogeneousBwdTrans(physTot, previous_iter[i],
625 wk[i]);
626 }
627 }
628 else
629 {
630 for (size_t i = 0; i < nvel; ++i)
631 {
632 Vmath::Vcopy(physTot, previous_iter[i], 1, wk[i], 1);
633 }
634 }
635
636 // (L(U^i) - 1/alpha*U^i_jj)
637 m_mapping->VelocityLaplacian(wk, F_corrected,
638 1.0 / m_viscousRelaxation);
639
640 if (m_fields[0]->GetWaveSpace())
641 {
642 for (size_t i = 0; i < nvel; ++i)
643 {
644 m_fields[0]->HomogeneousFwdTrans(physTot, F_corrected[i],
645 F_corrected[i]);
646 }
647 }
648 else
649 {
650 for (size_t i = 0; i < nvel; ++i)
651 {
652 Vmath::Vcopy(physTot, F_corrected[i], 1, F_corrected[i], 1);
653 }
654 }
655
656 // Loop velocity components
657 for (size_t i = 0; i < nvel; ++i)
658 {
659 // (-alpha*L(U^i) + U^i_jj)
660 Vmath::Smul(physTot, -1.0 * m_viscousRelaxation, F_corrected[i],
661 1, F_corrected[i], 1);
662 // F_corrected = alpha*F + (-alpha*L(U^i) + U^i_jj)
663 Vmath::Smul(physTot, m_viscousRelaxation, Forcing[i], 1, wk[0],
664 1);
665 Vmath::Vadd(physTot, wk[0], 1, F_corrected[i], 1,
666 F_corrected[i], 1);
667
668 //
669 // Solve System
670 //
671 m_fields[i]->HelmSolve(F_corrected[i],
672 m_fields[i]->UpdateCoeffs(), factors);
673 m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
674
675 //
676 // Test convergence
677 //
678 error = m_fields[i]->L2(outarray[i], previous_iter[i]);
679
680 if (forcing_L2[i] != 0)
681 {
682 if ((error / forcing_L2[i] >= m_viscousTolerance))
683 {
684 converged = false;
685 }
686 }
687 else
688 {
689 if (error >= m_viscousTolerance)
690 {
691 converged = false;
692 }
693 }
694 if (error > max_error)
695 {
696 max_error = error;
697 }
698
699 // Copy field to previous_iter
700 Vmath::Vcopy(physTot, outarray[i], 1, previous_iter[i], 1);
701 }
702 }
703 if (m_verbose && m_session->GetComm()->GetRank() == 0)
704 {
705 std::cout << " Velocity system (mapping) converged in " << s
706 << " iterations with error = " << max_error << std::endl;
707 }
708 }
709}
710
711/**
712 * Explicit terms of the mapping
713 */
715 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
717{
718 size_t physTot = m_fields[0]->GetTotPoints();
723 for (int i = 0; i < m_nConvectiveFields; ++i)
724 {
725 velPhys[i] = Array<OneD, NekDouble>(physTot, 0.0);
726 Forcing[i] = Array<OneD, NekDouble>(physTot, 0.0);
727 tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
728 }
729
730 // Get fields and store velocity in wavespace and physical space
731 if (m_fields[0]->GetWaveSpace())
732 {
733 for (int i = 0; i < m_nConvectiveFields; ++i)
734 {
735 vel[i] = inarray[i];
736 m_fields[0]->HomogeneousBwdTrans(physTot, vel[i], velPhys[i]);
737 }
738 }
739 else
740 {
741 for (int i = 0; i < m_nConvectiveFields; ++i)
742 {
743 vel[i] = inarray[i];
744 Vmath::Vcopy(physTot, inarray[i], 1, velPhys[i], 1);
745 }
746 }
747
748 // Advection contribution
750
751 // Time-derivative contribution
752 if (m_mapping->IsTimeDependent())
753 {
754 MappingAccelerationCorrection(vel, velPhys, tmp);
755 for (int i = 0; i < m_nConvectiveFields; ++i)
756 {
757 Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
758 }
759 }
760
761 // Pressure contribution
763 {
765 for (int i = 0; i < m_nConvectiveFields; ++i)
766 {
767 Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
768 }
769 }
770 // Viscous contribution
772 {
773 MappingViscousCorrection(velPhys, tmp);
774 for (int i = 0; i < m_nConvectiveFields; ++i)
775 {
776 Vmath::Smul(physTot, m_kinvis, tmp[i], 1, tmp[i], 1);
777 Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
778 }
779 }
780
781 // If necessary, transform to wavespace
782 if (m_fields[0]->GetWaveSpace())
783 {
784 for (int i = 0; i < m_nConvectiveFields; ++i)
785 {
786 m_fields[0]->HomogeneousFwdTrans(physTot, Forcing[i], Forcing[i]);
787 }
788 }
789
790 // Add to outarray
791 for (int i = 0; i < m_nConvectiveFields; ++i)
792 {
793 Vmath::Vadd(physTot, outarray[i], 1, Forcing[i], 1, outarray[i], 1);
794 }
795}
796
798 const Array<OneD, const Array<OneD, NekDouble>> &velPhys,
800{
801 size_t physTot = m_fields[0]->GetTotPoints();
802 size_t nvel = m_nConvectiveFields;
803
804 Array<OneD, Array<OneD, NekDouble>> wk(nvel * nvel);
805
806 // Apply Christoffel symbols to obtain {i,kj}vel(k)
807 m_mapping->ApplyChristoffelContravar(velPhys, wk);
808
809 // Calculate correction -U^j*{i,kj}vel(k)
810 for (size_t i = 0; i < nvel; i++)
811 {
812 Vmath::Zero(physTot, outarray[i], 1);
813 for (size_t j = 0; j < nvel; j++)
814 {
815 Vmath::Vvtvp(physTot, wk[i * nvel + j], 1, velPhys[j], 1,
816 outarray[i], 1, outarray[i], 1);
817 }
818 Vmath::Neg(physTot, outarray[i], 1);
819 }
820}
821
823 const Array<OneD, const Array<OneD, NekDouble>> &vel,
824 const Array<OneD, const Array<OneD, NekDouble>> &velPhys,
826{
827 size_t physTot = m_fields[0]->GetTotPoints();
828 size_t nvel = m_nConvectiveFields;
829
830 Array<OneD, Array<OneD, NekDouble>> wk(nvel * nvel);
833 for (size_t i = 0; i < nvel; i++)
834 {
835 tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
836 coordVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
837 }
838 // Get coordinates velocity in transformed system
839 m_mapping->GetCoordVelocity(tmp);
840 m_mapping->ContravarFromCartesian(tmp, coordVel);
841
842 // Calculate first term: U^j u^i,j = U^j (du^i/dx^j + {i,kj}u^k)
843 m_mapping->ApplyChristoffelContravar(velPhys, wk);
844 for (size_t i = 0; i < nvel; i++)
845 {
846 Vmath::Zero(physTot, outarray[i], 1);
847
848 m_fields[0]->PhysDeriv(velPhys[i], tmp[0], tmp[1]);
849 for (size_t j = 0; j < nvel; j++)
850 {
851 if (j == 2)
852 {
853 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j], vel[i],
854 tmp[2]);
855 if (m_fields[0]->GetWaveSpace())
856 {
857 m_fields[0]->HomogeneousBwdTrans(physTot, tmp[2], tmp[2]);
858 }
859 }
860
861 Vmath::Vadd(physTot, wk[i * nvel + j], 1, tmp[j], 1,
862 wk[i * nvel + j], 1);
863
864 Vmath::Vvtvp(physTot, coordVel[j], 1, wk[i * nvel + j], 1,
865 outarray[i], 1, outarray[i], 1);
866 }
867 }
868
869 // Set wavespace to false and store current value
870 bool wavespace = m_fields[0]->GetWaveSpace();
871 m_fields[0]->SetWaveSpace(false);
872
873 // Add -u^j U^i,j
874 m_mapping->ApplyChristoffelContravar(coordVel, wk);
875 for (size_t i = 0; i < nvel; i++)
876 {
877 if (nvel == 2)
878 {
879 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1]);
880 }
881 else
882 {
883 m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1], tmp[2]);
884 }
885
886 for (size_t j = 0; j < nvel; j++)
887 {
888 Vmath::Vadd(physTot, wk[i * nvel + j], 1, tmp[j], 1,
889 wk[i * nvel + j], 1);
890 Vmath::Neg(physTot, wk[i * nvel + j], 1);
891
892 Vmath::Vvtvp(physTot, velPhys[j], 1, wk[i * nvel + j], 1,
893 outarray[i], 1, outarray[i], 1);
894 }
895 }
896
897 // Restore value of wavespace
898 m_fields[0]->SetWaveSpace(wavespace);
899}
900
903{
904 size_t physTot = m_fields[0]->GetTotPoints();
905 size_t nvel = m_nConvectiveFields;
906
907 // Calculate g^(ij)p_(,j)
908 m_mapping->RaiseIndex(m_gradP, outarray);
909
910 // Calculate correction = (nabla p)/J - g^(ij)p_,j
911 // (Jac is not required if it is constant)
912 if (!m_mapping->HasConstantJacobian())
913 {
914 Array<OneD, NekDouble> Jac(physTot, 0.0);
915 m_mapping->GetJacobian(Jac);
916 for (size_t i = 0; i < nvel; ++i)
917 {
918 Vmath::Vdiv(physTot, m_gradP[i], 1, Jac, 1, m_gradP[i], 1);
919 }
920 }
921 for (size_t i = 0; i < nvel; ++i)
922 {
923 Vmath::Vsub(physTot, m_gradP[i], 1, outarray[i], 1, outarray[i], 1);
924 }
925}
926
928 const Array<OneD, const Array<OneD, NekDouble>> &velPhys,
930{
931 // L(U) - 1.0*d^2(u^i)/dx^jdx^j
932 m_mapping->VelocityLaplacian(velPhys, outarray, 1.0);
933}
934
935} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:264
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
NekDouble m_kinvis
Kinematic viscosity.
bool m_SmoothAdvection
bool to identify if advection term smoothing is requested
ExtrapolateSharedPtr m_extrapolation
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
int m_intSteps
Number of time integration steps AND Order of extrapolation for pressure boundary conditions.
int m_nConvectiveFields
Number of fields to be convected;.
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
NekDouble m_timestep
Time step size.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTotPoints()
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:71
Base class for unsteady solvers.
SOLVER_UTILS_EXPORT void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
LibUtilities::TimeIntegrationSchemeSharedPtr m_intScheme
Wrapper to the time integration scheme.
void v_SolvePressure(const Array< OneD, NekDouble > &Forcing) override
Definition: VCSMapping.cpp:412
void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt) override
Definition: VCSMapping.cpp:210
void MappingViscousCorrection(const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: VCSMapping.cpp:927
void v_DoInitialise(bool dumpInitialConditions=true) override
Sets up initial conditions.
Definition: VCSMapping.cpp:128
void v_SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt) override
Definition: VCSMapping.cpp:335
static std::string className
Name of class.
Definition: VCSMapping.h:61
NekDouble m_pressureTolerance
Definition: VCSMapping.h:82
NekDouble m_viscousRelaxation
Definition: VCSMapping.h:85
NekDouble m_pressureRelaxation
Definition: VCSMapping.h:84
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Definition: VCSMapping.h:50
NekDouble m_viscousTolerance
Definition: VCSMapping.h:83
void v_EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time) override
Definition: VCSMapping.cpp:167
void v_SolveViscous(const Array< OneD, const Array< OneD, NekDouble > > &Forcing, const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble aii_Dt) override
Definition: VCSMapping.cpp:550
GlobalMapping::MappingSharedPtr m_mapping
Definition: VCSMapping.h:69
void MappingPressureCorrection(Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: VCSMapping.cpp:901
void ApplyIncNSMappingForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: VCSMapping.cpp:714
Array< OneD, Array< OneD, NekDouble > > m_presForcingCorrection
Definition: VCSMapping.h:124
VCSMapping(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Definition: VCSMapping.cpp:58
static std::string solverTypeLookupId
Definition: VCSMapping.h:71
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
Definition: VCSMapping.cpp:65
void MappingAccelerationCorrection(const Array< OneD, const Array< OneD, NekDouble > > &vel, const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: VCSMapping.cpp:822
Array< OneD, Array< OneD, NekDouble > > m_gradP
Definition: VCSMapping.h:88
void MappingAdvectionCorrection(const Array< OneD, const Array< OneD, NekDouble > > &velPhys, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: VCSMapping.cpp:797
NekDouble m_sVVCutoffRatio
cutt off ratio from which to start decayhing modes
NekDouble m_sVVDiffCoeff
Diffusion coefficient of SVV modes.
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
Array< OneD, Array< OneD, NekDouble > > m_F
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
virtual void v_SolveViscous(const Array< OneD, const Array< OneD, NekDouble > > &Forcing, const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble aii_Dt)
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:149
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:135
std::shared_ptr< SessionReader > SessionReaderSharedPtr
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:87
EquationSystemFactory & GetEquationSystemFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
StdRegions::ConstFactorMap factors
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:47
double NekDouble
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 Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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
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 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