Nektar++
CompressibleFlowSystemImplicit.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File CFSImplicit.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: Compressible flow system base class with auxiliary functions
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
38 
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 CFSImplicit::CFSImplicit(const LibUtilities::SessionReaderSharedPtr &pSession,
47  : UnsteadySystem(pSession, pGraph), CompressibleFlowSystem(pSession, pGraph)
48 {
49 }
50 
51 /**
52  * @brief Initialization object for CFSImplicit class.
53  */
54 void CFSImplicit::v_InitObject(bool DeclareFields)
55 {
57  m_explicitAdvection = false;
58  m_explicitDiffusion = false;
59 
60  // initialise implicit parameters
61  m_session->LoadParameter("JacobiFreeEps", m_jacobiFreeEps, 5.0E-8);
62 
63  int ntmp;
64  m_session->LoadParameter("AdvectionJacFlag", ntmp, 1);
65  m_advectionJacFlag = (0 == ntmp) ? false : true;
66 
67  m_session->LoadParameter("ViscousJacFlag", ntmp, 1);
68  m_viscousJacFlag = (0 == ntmp) ? false : true;
69 
70  // initialise implicit functors
72 
74 
75  m_flagImplicitSolver = true;
76 }
77 
79 {
80  std::string SolverType = "Newton";
81  if (m_session->DefinesSolverInfo("NonlinSysIterSolver"))
82  {
83  SolverType = m_session->GetSolverInfo("NonlinSysIterSolver");
84  }
85  ASSERTL0(LibUtilities::GetNekNonlinSysFactory().ModuleExists(SolverType),
86  "NekNonlinSys '" + SolverType + "' is not defined.\n");
87  int ntotal = m_fields[0]->GetNcoeffs() * m_fields.size();
88 
90 
91  key.m_NonlinIterTolRelativeL2 = 1.0E-3;
92  key.m_LinSysRelativeTolInNonlin = 5.0E-2;
94  key.m_NekLinSysMaxIterations = 30;
95  key.m_LinSysMaxStorage = 30;
96 
98  SolverType, m_session, m_comm, ntotal, key);
99 
103  this);
105 
106  int nvariables = m_fields.size();
107  const MultiRegions::LocTraceToTraceMapSharedPtr locTraceToTraceMap =
108  m_fields[0]->GetLocTraceToTraceMap();
109 
110  locTraceToTraceMap->CalcLocTracePhysToTraceIDMap(m_fields[0]->GetTrace(),
111  m_spacedim);
112  for (int i = 1; i < nvariables; i++)
113  {
114  m_fields[i]->GetLocTraceToTraceMap()->SetLocTracePhysToTraceIDMap(
115  locTraceToTraceMap->GetLocTracephysToTraceIDMap());
116  }
117 
119  "PreconCfsBRJ", m_fields, m_session, m_comm);
120  NekPreconCfsOperators tmpPreconOp;
122  this);
123  m_preconCfs->SetOperators(tmpPreconOp);
124 
125  m_session->LoadParameter("NewtonAbsoluteIteTol", m_newtonAbsoluteIteTol,
126  1.0E-12);
127 
128  m_session->LoadParameter("nPadding", m_nPadding, 4);
129 
130  m_nonlinsol->SetSysOperators(nekSysOp);
131 }
132 
133 /**
134  * @brief Destructor for CFSImplicit class.
135  */
137 {
138 }
139 
142  const bool &flag)
143 {
144  const Array<OneD, const NekDouble> refsource =
145  m_nonlinsol->GetRefSourceVec();
146  NonlinSysEvaluatorCoeff(inarray, out, flag, refsource);
147 }
148 
151  const bool &flag, const Array<OneD, const NekDouble> &source)
152 {
153  LibUtilities::Timer timer;
154  boost::ignore_unused(flag);
155  unsigned int nvariables = m_fields.size();
156  unsigned int npoints = m_fields[0]->GetNcoeffs();
157  Array<OneD, Array<OneD, NekDouble>> in2D(nvariables);
158  Array<OneD, Array<OneD, NekDouble>> out2D(nvariables);
159  Array<OneD, Array<OneD, NekDouble>> source2D(nvariables);
160  for (int i = 0; i < nvariables; ++i)
161  {
162  int offset = i * npoints;
163  in2D[i] = inarray + offset;
164  out2D[i] = out + offset;
165  source2D[i] = source + offset;
166  }
167 
168  timer.Start();
169  NonlinSysEvaluatorCoeff(in2D, out2D, source2D);
170  timer.Stop();
171  timer.AccumulateRegion("CFSImplicit::NonlinSysEvaluatorCoeff");
172 }
173 
175  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
177  const Array<OneD, const Array<OneD, NekDouble>> &source)
178 {
179  LibUtilities::Timer timer;
180  unsigned int nvariable = inarray.size();
181  unsigned int ncoeffs = inarray[nvariable - 1].size();
182  unsigned int npoints = m_fields[0]->GetNpoints();
183 
184  Array<OneD, Array<OneD, NekDouble>> inpnts(nvariable);
185 
186  for (int i = 0; i < nvariable; ++i)
187  {
188  inpnts[i] = Array<OneD, NekDouble>(npoints, 0.0);
189  m_fields[i]->BwdTrans(inarray[i], inpnts[i]);
190  }
191 
192  timer.Start();
193  DoOdeProjection(inpnts, inpnts, m_bndEvaluateTime);
194  timer.Stop();
195  timer.AccumulateRegion("CompressibleFlowSystem::DoOdeProjection", 1);
196 
197  timer.Start();
198  DoOdeRhsCoeff(inpnts, out, m_bndEvaluateTime);
199  timer.Stop();
200  timer.AccumulateRegion("CFSImplicit::DoOdeRhsCoeff", 1);
201 
202  for (int i = 0; i < nvariable; ++i)
203  {
204  Vmath::Svtvp(ncoeffs, -m_TimeIntegLambda, out[i], 1, inarray[i], 1,
205  out[i], 1);
206  }
207 
208  if (NullNekDoubleArrayOfArray != source)
209  {
210  for (int i = 0; i < nvariable; ++i)
211  {
212  Vmath::Vsub(ncoeffs, out[i], 1, source[i], 1, out[i], 1);
213  }
214  }
215  return;
216 }
217 
218 /**
219  * @brief Compute the right-hand side.
220  */
222  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
223  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
224 {
225  LibUtilities::Timer timer;
226 
227  int nvariables = inarray.size();
228  int nTracePts = GetTraceTotPoints();
229  int ncoeffs = GetNcoeffs();
230  // Store forwards/backwards space along trace space
231  Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
232  Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
233 
235  {
238  }
239  else
240  {
241  for (int i = 0; i < nvariables; ++i)
242  {
243  Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
244  Bwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
245  m_fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
246  }
247  }
248 
249  // Calculate advection
250  timer.Start();
251  DoAdvectionCoeff(inarray, outarray, time, Fwd, Bwd);
252  timer.Stop();
253  timer.AccumulateRegion("CFSImplicit::DoAdvectionCoeff", 2);
254 
255  // Negate results
256  for (int i = 0; i < nvariables; ++i)
257  {
258  Vmath::Neg(ncoeffs, outarray[i], 1);
259  }
260 
261  timer.Start();
262  DoDiffusionCoeff(inarray, outarray, Fwd, Bwd);
263  timer.Stop();
264  timer.AccumulateRegion("CFSImplicit::DoDiffusionCoeff", 2);
265 
266  // Add forcing terms
267  for (auto &x : m_forcing)
268  {
269  x->ApplyCoeff(m_fields, inarray, outarray, time);
270  }
271 
272  if (m_useLocalTimeStep)
273  {
274  int nElements = m_fields[0]->GetExpSize();
275  int nq, offset;
276  NekDouble fac;
278 
279  Array<OneD, NekDouble> tstep(nElements, 0.0);
280  GetElmtTimeStep(inarray, tstep);
281 
282  // Loop over elements
283  for (int n = 0; n < nElements; ++n)
284  {
285  nq = m_fields[0]->GetExp(n)->GetNcoeffs();
286  offset = m_fields[0]->GetCoeff_Offset(n);
287  fac = tstep[n] / m_timestep;
288  for (int i = 0; i < nvariables; ++i)
289  {
290  Vmath::Smul(nq, fac, outarray[i] + offset, 1,
291  tmp = outarray[i] + offset, 1);
292  }
293  }
294  }
295 }
296 
297 /**
298  * @brief Compute the advection terms for the right-hand side
299  */
301  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
302  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time,
303  const Array<OneD, const Array<OneD, NekDouble>> &pFwd,
304  const Array<OneD, const Array<OneD, NekDouble>> &pBwd)
305 {
306  int nvariables = inarray.size();
308 
309  m_advObject->AdvectCoeffs(nvariables, m_fields, advVel, inarray, outarray,
310  time, pFwd, pBwd);
311 }
312 
314  const Array<OneD, const Array<OneD, NekDouble>> &inpnts,
315  Array<OneD, Array<OneD, NekDouble>> &outpnt, const NekDouble time,
316  const NekDouble lambda)
317 {
318  unsigned int nvariables = inpnts.size();
319  unsigned int ncoeffs = m_fields[0]->GetNcoeffs();
320  unsigned int ntotal = nvariables * ncoeffs;
321 
322  Array<OneD, NekDouble> inarray(ntotal);
323  Array<OneD, NekDouble> out(ntotal);
324  Array<OneD, NekDouble> tmpArray;
325 
326  for (int i = 0; i < nvariables; ++i)
327  {
328  int noffset = i * ncoeffs;
329  tmpArray = inarray + noffset;
330  m_fields[i]->FwdTrans(inpnts[i], tmpArray);
331  }
332 
333  DoImplicitSolveCoeff(inpnts, inarray, out, time, lambda);
334 
335  for (int i = 0; i < nvariables; ++i)
336  {
337  int noffset = i * ncoeffs;
338  tmpArray = out + noffset;
339  m_fields[i]->BwdTrans(tmpArray, outpnt[i]);
340  }
341 }
342 
344  const Array<OneD, const Array<OneD, NekDouble>> &inpnts,
346  const NekDouble time, const NekDouble lambda)
347 {
348  boost::ignore_unused(inpnts);
349 
350  m_TimeIntegLambda = lambda;
351  m_bndEvaluateTime = time;
352  m_solutionPhys = inpnts;
353  unsigned int ntotal = inarray.size();
354 
355  if (m_inArrayNorm < 0.0)
356  {
357  CalcRefValues(inarray);
358  }
359 
360  NekDouble tol2 =
362 
363  m_nonlinsol->v_SetupNekNonlinSystem(ntotal, inarray, inarray, 0);
364 
365  m_TotNewtonIts += m_nonlinsol->SolveSystem(ntotal, inarray, out, 0, tol2);
366 
367  m_TotLinIts += m_nonlinsol->GetNtotLinSysIts();
368 
369  m_TotImpStages++;
370  m_StagesPerStep++;
371 }
372 
374 {
375  unsigned int nvariables = m_fields.size();
376  unsigned int ntotal = inarray.size();
377  unsigned int npoints = ntotal / nvariables;
378 
379  unsigned int nTotalGlobal = ntotal;
380  m_comm->AllReduce(nTotalGlobal, Nektar::LibUtilities::ReduceSum);
381  unsigned int nTotalDOF = nTotalGlobal / nvariables;
382  NekDouble invTotalDOF = 1.0 / nTotalDOF;
383 
384  m_inArrayNorm = 0.0;
385  m_magnitdEstimat = Array<OneD, NekDouble>(nvariables, 0.0);
386 
387  for (int i = 0; i < nvariables; ++i)
388  {
389  int offset = i * npoints;
390  m_magnitdEstimat[i] =
391  Vmath::Dot(npoints, inarray + offset, inarray + offset);
392  }
394 
395  for (int i = 0; i < nvariables; ++i)
396  {
398  }
399 
400  for (int i = 2; i < nvariables - 1; ++i)
401  {
403  }
404  for (int i = 2; i < nvariables - 1; ++i)
405  {
407  }
408 
409  for (int i = 0; i < nvariables; ++i)
410  {
411  m_magnitdEstimat[i] = sqrt(m_magnitdEstimat[i] * invTotalDOF);
412  }
413  if (m_root && m_verbose)
414  {
415  for (int i = 0; i < nvariables; ++i)
416  {
417  cout << "m_magnitdEstimat[" << i << "] = " << m_magnitdEstimat[i]
418  << endl;
419  }
420  cout << "m_inArrayNorm = " << m_inArrayNorm << endl;
421  }
422 }
423 
425  Array<OneD, NekDouble> &outarray,
426  const bool &flag)
427 {
428  LibUtilities::Timer timer, Gtimer;
429 
430  Gtimer.Start();
431  if (m_preconCfs->UpdatePreconMatCheck(NullNekDouble1DArray,
434  {
435  int nvariables = m_solutionPhys.size();
436 
437  int nphspnt = m_solutionPhys[nvariables - 1].size();
438  Array<OneD, Array<OneD, NekDouble>> intmp(nvariables);
439  for (int i = 0; i < nvariables; i++)
440  {
441  intmp[i] = Array<OneD, NekDouble>(nphspnt, 0.0);
442  }
443 
444  timer.Start();
446  timer.Stop();
447  timer.AccumulateRegion("CompressibleFlowSystem::DoOdeProjection", 1);
448 
449  timer.Start();
450  m_preconCfs->BuildPreconCfs(m_fields, intmp, m_bndEvaluateTime,
452  timer.Stop();
453  timer.AccumulateRegion("PreconCfsOp::BuildPreconCfs", 1);
454  }
455 
456  m_flagUpdatePreconMat = false;
457 
458  timer.Start();
459  m_preconCfs->DoPreconCfs(m_fields, inarray, outarray, flag);
460  timer.Stop();
461  timer.AccumulateRegion("PreconCfsOp::DoPreconCfs", 1);
462 
463  Gtimer.Stop();
464  Gtimer.AccumulateRegion("CFSImplicit::PreconCoeff");
465 }
466 
467 template <typename DataType, typename TypeNekBlkMatSharedPtr>
469  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
470  const Array<OneD, const TensorOfArray2D<NekDouble>> &qfield,
472  TensorOfArray4D<DataType> &StdMatDataDBB,
473  TensorOfArray5D<DataType> &StdMatDataDBDB)
474 {
475  if (StdMatDataDBB.size() == 0)
476  {
477  CalcVolJacStdMat(StdMatDataDBB, StdMatDataDBDB);
478  }
479 
480  int nSpaceDim = m_graph->GetSpaceDimension();
481  int nvariable = inarray.size();
482  int npoints = m_fields[0]->GetTotPoints();
483  int nVar2 = nvariable * nvariable;
484  std::shared_ptr<LocalRegions::ExpansionVector> expvect =
485  m_fields[0]->GetExp();
486  int nTotElmt = (*expvect).size();
487 
488  Array<OneD, NekDouble> mu(npoints, 0.0);
489  Array<OneD, NekDouble> DmuDT(npoints, 0.0);
490  if (m_viscousJacFlag)
491  {
492  CalcMuDmuDT(inarray, mu, DmuDT);
493  }
494 
495  Array<OneD, NekDouble> normals;
497  for (int i = 0; i < 3; i++)
498  {
499  normal3D[i] = Array<OneD, NekDouble>(3, 0.0);
500  }
501  normal3D[0][0] = 1.0;
502  normal3D[1][1] = 1.0;
503  normal3D[2][2] = 1.0;
505 
506  DNekMatSharedPtr wspMat =
507  MemoryManager<DNekMat>::AllocateSharedPtr(nvariable, nvariable, 0.0);
509  nvariable - 1, nvariable, 0.0);
510 
511  Array<OneD, DataType> GmatxData;
512  Array<OneD, DataType> MatData;
513 
514  Array<OneD, NekDouble> tmppnts;
515  TensorOfArray3D<NekDouble> PntJacCons(
516  m_spacedim); // Nvar*Nvar*Ndir*Nelmt*Npnt
517  TensorOfArray3D<DataType> PntJacConsStd(
518  m_spacedim); // Nvar*Nvar*Ndir*Nelmt*Npnt
521  TensorOfArray4D<NekDouble> PntJacDerv(
522  m_spacedim); // Nvar*Nvar*Ndir*Nelmt*Npnt
523  TensorOfArray4D<DataType> PntJacDervStd(
524  m_spacedim); // Nvar*Nvar*Ndir*Nelmt*Npnt
526  m_spacedim); // Nvar*Nvar*Ndir*Nelmt*Npnt
528  m_spacedim); // Nvar*Nvar*Ndir*Nelmt*Npnt
529  for (int ndir = 0; ndir < m_spacedim; ndir++)
530  {
531  PntJacDerv[ndir] = TensorOfArray3D<NekDouble>(m_spacedim);
532  PntJacDervStd[ndir] = TensorOfArray3D<DataType>(m_spacedim);
535  }
536 
538  Array<OneD, NekDouble> locDmuDT;
539  Array<OneD, Array<OneD, NekDouble>> locVars(nvariable);
541  for (int ndir = 0; ndir < m_spacedim; ndir++)
542  {
543  locDerv[ndir] = Array<OneD, Array<OneD, NekDouble>>(nvariable);
544  }
545 
546  int nElmtCoefOld = -1;
547  for (int ne = 0; ne < nTotElmt; ne++)
548  {
549  int nElmtCoef = (*expvect)[ne]->GetNcoeffs();
550  int nElmtCoef2 = nElmtCoef * nElmtCoef;
551  int nElmtPnt = (*expvect)[ne]->GetTotPoints();
552 
553  int nQuot = nElmtCoef2 / m_nPadding;
554  int nRemd = nElmtCoef2 - nQuot * m_nPadding;
555  int nQuotPlus = nQuot;
556  if (nRemd > 0)
557  {
558  nQuotPlus++;
559  }
560  int nElmtCoef2Paded = nQuotPlus * m_nPadding;
561 
562  if (nElmtPnt > PntJacCons[0].size() || nElmtCoef > nElmtCoefOld)
563  {
564  nElmtCoefOld = nElmtCoef;
565  for (int ndir = 0; ndir < 3; ndir++)
566  {
567  normalPnt[ndir] = Array<OneD, NekDouble>(npoints, 0.0);
568  }
569  tmppnts = Array<OneD, NekDouble>(nElmtPnt);
570  MatData = Array<OneD, DataType>(nElmtCoef2Paded * nVar2);
571  for (int ndir = 0; ndir < m_spacedim; ndir++)
572  {
573  ConsCurv[ndir] = Array<OneD, NekDouble>(nElmtPnt);
574  ConsStdd[ndir] = Array<OneD, NekDouble>(nElmtPnt);
575  PntJacCons[ndir] =
577  PntJacConsStd[ndir] =
579  for (int i = 0; i < nElmtPnt; i++)
580  {
581  PntJacCons[ndir][i] = Array<OneD, NekDouble>(nVar2);
582  PntJacConsStd[ndir][i] = Array<OneD, DataType>(nVar2);
583  }
584 
585  for (int ndir1 = 0; ndir1 < m_spacedim; ndir1++)
586  {
587  PntJacDerv[ndir][ndir1] =
589  PntJacDervStd[ndir][ndir1] =
591  DervStdd[ndir][ndir1] = Array<OneD, NekDouble>(nElmtPnt);
592  DervCurv[ndir][ndir1] = Array<OneD, NekDouble>(nElmtPnt);
593  for (int i = 0; i < nElmtPnt; i++)
594  {
595  PntJacDerv[ndir][ndir1][i] =
596  Array<OneD, NekDouble>(nVar2);
597  PntJacDervStd[ndir][ndir1][i] =
598  Array<OneD, DataType>(nVar2);
599  }
600  }
601  }
602  }
603 
604  int noffset = GetPhys_Offset(ne);
605  for (int j = 0; j < nvariable; j++)
606  {
607  locVars[j] = inarray[j] + noffset;
608  }
609 
610  if (m_advectionJacFlag)
611  {
612  for (int nFluxDir = 0; nFluxDir < nSpaceDim; nFluxDir++)
613  {
614  normals = normal3D[nFluxDir];
615  GetFluxVectorJacDirElmt(nvariable, nElmtPnt, locVars, normals,
616  wspMat, PntJacCons[nFluxDir]);
617  }
618  }
619 
620  if (m_viscousJacFlag)
621  {
622  for (int j = 0; j < nSpaceDim; j++)
623  {
624  for (int k = 0; k < nvariable; k++)
625  {
626  locDerv[j][k] = qfield[j][k] + noffset;
627  }
628  }
629  locmu = mu + noffset;
630  locDmuDT = DmuDT + noffset;
631  for (int nFluxDir = 0; nFluxDir < nSpaceDim; nFluxDir++)
632  {
633  normals = normal3D[nFluxDir];
634  MinusDiffusionFluxJacPoint(nvariable, nElmtPnt, locVars,
635  locDerv, locmu, locDmuDT, normals,
636  wspMatDrv, PntJacCons[nFluxDir]);
637  }
638  }
639 
640  if (m_viscousJacFlag)
641  {
642  locmu = mu + noffset;
643  for (int nFluxDir = 0; nFluxDir < nSpaceDim; nFluxDir++)
644  {
645  Vmath::Fill(npoints, 1.0, normalPnt[nFluxDir], 1);
646  for (int nDervDir = 0; nDervDir < nSpaceDim; nDervDir++)
647  {
649  nvariable, nElmtPnt, nDervDir, locVars, locmu,
650  normalPnt, wspMatDrv, PntJacDerv[nFluxDir][nDervDir]);
651  }
652  Vmath::Fill(npoints, 0.0, normalPnt[nFluxDir], 1);
653  }
654  }
655 
656  for (int n = 0; n < nvariable; n++)
657  {
658  for (int m = 0; m < nvariable; m++)
659  {
660  int nVarOffset = m + n * nvariable;
661  GmatxData = gmtxarray[m][n]->GetBlock(ne, ne)->GetPtr();
662 
663  for (int ndStd0 = 0; ndStd0 < m_spacedim; ndStd0++)
664  {
665  Vmath::Zero(nElmtPnt, ConsStdd[ndStd0], 1);
666  }
667  for (int ndir = 0; ndir < m_spacedim; ndir++)
668  {
669  for (int i = 0; i < nElmtPnt; i++)
670  {
671  tmppnts[i] = PntJacCons[ndir][i][nVarOffset];
672  }
673  (*expvect)[ne]->AlignVectorToCollapsedDir(ndir, tmppnts,
674  ConsCurv);
675  for (int nd = 0; nd < m_spacedim; nd++)
676  {
677  Vmath::Vadd(nElmtPnt, ConsCurv[nd], 1, ConsStdd[nd], 1,
678  ConsStdd[nd], 1);
679  }
680  }
681 
682  for (int ndir = 0; ndir < m_spacedim; ndir++)
683  {
684  (*expvect)[ne]->MultiplyByQuadratureMetric(
685  ConsStdd[ndir], ConsStdd[ndir]); // weight with metric
686  for (int i = 0; i < nElmtPnt; i++)
687  {
688  PntJacConsStd[ndir][i][nVarOffset] =
689  DataType(ConsStdd[ndir][i]);
690  }
691  }
692  }
693  }
694 
695  if (m_viscousJacFlag)
696  {
697  for (int m = 0; m < nvariable; m++)
698  {
699  for (int n = 0; n < nvariable; n++)
700  {
701  int nVarOffset = m + n * nvariable;
702  for (int ndStd0 = 0; ndStd0 < m_spacedim; ndStd0++)
703  {
704  for (int ndStd1 = 0; ndStd1 < m_spacedim; ndStd1++)
705  {
706  Vmath::Zero(nElmtPnt, DervStdd[ndStd0][ndStd1], 1);
707  }
708  }
709  for (int nd0 = 0; nd0 < m_spacedim; nd0++)
710  {
711  for (int nd1 = 0; nd1 < m_spacedim; nd1++)
712  {
713  for (int i = 0; i < nElmtPnt; i++)
714  {
715  tmppnts[i] =
716  PntJacDerv[nd0][nd1][i][nVarOffset];
717  }
718 
719  (*expvect)[ne]->AlignVectorToCollapsedDir(
720  nd0, tmppnts, ConsCurv);
721  for (int nd = 0; nd < m_spacedim; nd++)
722  {
723  (*expvect)[ne]->AlignVectorToCollapsedDir(
724  nd1, ConsCurv[nd], DervCurv[nd]);
725  }
726 
727  for (int ndStd0 = 0; ndStd0 < m_spacedim; ndStd0++)
728  {
729  for (int ndStd1 = 0; ndStd1 < m_spacedim;
730  ndStd1++)
731  {
732  Vmath::Vadd(nElmtPnt,
733  DervCurv[ndStd0][ndStd1], 1,
734  DervStdd[ndStd0][ndStd1], 1,
735  DervStdd[ndStd0][ndStd1], 1);
736  }
737  }
738  }
739  }
740  for (int nd0 = 0; nd0 < m_spacedim; nd0++)
741  {
742  for (int nd1 = 0; nd1 < m_spacedim; nd1++)
743  {
744  (*expvect)[ne]->MultiplyByQuadratureMetric(
745  DervStdd[nd0][nd1],
746  DervStdd[nd0][nd1]); // weight with metric
747  for (int i = 0; i < nElmtPnt; i++)
748  {
749  PntJacDervStd[nd0][nd1][i][nVarOffset] =
750  -DataType(DervStdd[nd0][nd1][i]);
751  }
752  }
753  }
754  }
755  }
756  }
757 
758  Vmath::Zero(nElmtCoef2Paded * nVar2, MatData, 1);
759  DataType one = 1.0;
760  for (int ndir = 0; ndir < m_spacedim; ndir++)
761  {
762  for (int i = 0; i < nElmtPnt; i++)
763  {
764  Blas::Ger(nElmtCoef2Paded, nVar2, one,
765  &StdMatDataDBB[ne][ndir][i][0], 1,
766  &PntJacConsStd[ndir][i][0], 1, &MatData[0],
767  nElmtCoef2Paded);
768  }
769  }
770 
771  if (m_viscousJacFlag)
772  {
773  for (int nd0 = 0; nd0 < m_spacedim; nd0++)
774  {
775  for (int nd1 = 0; nd1 < m_spacedim; nd1++)
776  {
777  for (int i = 0; i < nElmtPnt; i++)
778  {
779  Blas::Ger(nElmtCoef2Paded, nVar2, one,
780  &StdMatDataDBDB[ne][nd0][nd1][i][0], 1,
781  &PntJacDervStd[nd0][nd1][i][0], 1,
782  &MatData[0], nElmtCoef2Paded);
783  }
784  }
785  }
786  }
787 
789 
790  for (int n = 0; n < nvariable; n++)
791  {
792  for (int m = 0; m < nvariable; m++)
793  {
794  int nVarOffset = m + n * nvariable;
795  GmatxData = gmtxarray[m][n]->GetBlock(ne, ne)->GetPtr();
796  Vmath::Vcopy(nElmtCoef2,
797  tmpA = MatData + nVarOffset * nElmtCoef2Paded, 1,
798  GmatxData, 1);
799  }
800  }
801  }
802 }
803 
804 template <typename DataType>
806  TensorOfArray5D<DataType> &StdMatDataDBDB)
807 {
808  std::shared_ptr<LocalRegions::ExpansionVector> expvect =
809  m_fields[0]->GetExp();
810  int nTotElmt = (*expvect).size();
811 
812  StdMatDataDBB = TensorOfArray4D<DataType>(nTotElmt);
813  StdMatDataDBDB = TensorOfArray5D<DataType>(nTotElmt);
814 
815  vector<DNekMatSharedPtr> VectStdDerivBase0;
816  vector<TensorOfArray3D<DataType>> VectStdDerivBase_Base;
817  vector<TensorOfArray4D<DataType>> VectStdDervBase_DervBase;
818  DNekMatSharedPtr MatStdDerivBase0;
821  for (int ne = 0; ne < nTotElmt; ne++)
822  {
824  stdExp = (*expvect)[ne]->GetStdExp();
826  stdExp->DetShapeType(), *stdExp);
827  MatStdDerivBase0 = stdExp->GetStdMatrix(matkey);
828 
829  int nTotStdExp = VectStdDerivBase0.size();
830  int nFoundStdExp = -1;
831  for (int i = 0; i < nTotStdExp; i++)
832  {
833  if ((*VectStdDerivBase0[i]) == (*MatStdDerivBase0))
834  {
835  nFoundStdExp = i;
836  }
837  }
838  if (nFoundStdExp >= 0)
839  {
840  StdMatDataDBB[ne] = VectStdDerivBase_Base[nFoundStdExp];
841  StdMatDataDBDB[ne] = VectStdDervBase_DervBase[nFoundStdExp];
842  }
843  else
844  {
845  int nElmtCoef = (*expvect)[ne]->GetNcoeffs();
846  int nElmtCoef2 = nElmtCoef * nElmtCoef;
847  int nElmtPnt = (*expvect)[ne]->GetTotPoints();
848 
849  int nQuot = nElmtCoef2 / m_nPadding;
850  int nRemd = nElmtCoef2 - nQuot * m_nPadding;
851  int nQuotPlus = nQuot;
852  if (nRemd > 0)
853  {
854  nQuotPlus++;
855  }
856  int nPaded = nQuotPlus * m_nPadding;
857 
858  ArrayStdMat[0] = MatStdDerivBase0;
859  if (m_spacedim > 1)
860  {
862  StdRegions::eDerivBase1, stdExp->DetShapeType(), *stdExp);
863  ArrayStdMat[1] = stdExp->GetStdMatrix(matkey);
864 
865  if (m_spacedim > 2)
866  {
868  stdExp->DetShapeType(),
869  *stdExp);
870  ArrayStdMat[2] = stdExp->GetStdMatrix(matkey);
871  }
872  }
873  for (int nd0 = 0; nd0 < m_spacedim; nd0++)
874  {
875  ArrayStdMatData[nd0] = ArrayStdMat[nd0]->GetPtr();
876  }
877 
879  stdExp->DetShapeType(), *stdExp);
880  DNekMatSharedPtr BwdMat = stdExp->GetStdMatrix(matkey);
881  Array<OneD, NekDouble> BwdMatData = BwdMat->GetPtr();
882 
885 
886  for (int nd0 = 0; nd0 < m_spacedim; nd0++)
887  {
888  tmpStdDBB[nd0] = Array<OneD, Array<OneD, DataType>>(nElmtPnt);
889  for (int i = 0; i < nElmtPnt; i++)
890  {
891  tmpStdDBB[nd0][i] = Array<OneD, DataType>(nPaded, 0.0);
892  for (int nc1 = 0; nc1 < nElmtCoef; nc1++)
893  {
894  int noffset = nc1 * nElmtCoef;
895  for (int nc0 = 0; nc0 < nElmtCoef; nc0++)
896  {
897  tmpStdDBB[nd0][i][nc0 + noffset] = DataType(
898  ArrayStdMatData[nd0][i * nElmtCoef + nc0] *
899  BwdMatData[i * nElmtCoef + nc1]);
900  }
901  }
902  }
903 
904  tmpStdDBDB[nd0] = TensorOfArray3D<DataType>(m_spacedim);
905  for (int nd1 = 0; nd1 < m_spacedim; nd1++)
906  {
907  tmpStdDBDB[nd0][nd1] =
909  for (int i = 0; i < nElmtPnt; i++)
910  {
911  tmpStdDBDB[nd0][nd1][i] =
912  Array<OneD, DataType>(nPaded, 0.0);
913  for (int nc1 = 0; nc1 < nElmtCoef; nc1++)
914  {
915  int noffset = nc1 * nElmtCoef;
916  for (int nc0 = 0; nc0 < nElmtCoef; nc0++)
917  {
918  tmpStdDBDB[nd0][nd1][i][nc0 + noffset] =
919  DataType(
920  ArrayStdMatData[nd0]
921  [i * nElmtCoef + nc0] *
922  ArrayStdMatData[nd1]
923  [i * nElmtCoef + nc1]);
924  }
925  }
926  }
927  }
928  }
929  VectStdDerivBase0.push_back(MatStdDerivBase0);
930  VectStdDerivBase_Base.push_back(tmpStdDBB);
931  VectStdDervBase_DervBase.push_back(tmpStdDBDB);
932 
933  StdMatDataDBB[ne] = tmpStdDBB;
934  StdMatDataDBDB[ne] = tmpStdDBDB;
935  }
936  }
937 }
938 
939 template <typename DataType, typename TypeNekBlkMatSharedPtr>
941  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
946  Array<OneD, Array<OneD, DataType>> &TraceJacDerivSign,
947  TensorOfArray5D<DataType> &TraceIPSymJacArray)
948 {
949  int nvariables = inarray.size();
950 
951  LibUtilities::Timer timer;
952  timer.Start();
953  GetTraceJac(inarray, qfield, TraceJac, TraceJacDeriv, TraceJacDerivSign,
954  TraceIPSymJacArray);
955  timer.Stop();
956  timer.AccumulateRegion("CFSImplicit::GetTraceJac", 10);
957 
960 
961  timer.Start();
962  m_advObject->AddTraceJacToMat(nvariables, m_spacedim, m_fields, TraceJac,
963  gmtxarray, tmpJac, tmpSign);
964  timer.Stop();
965  timer.AccumulateRegion("Advection::AddTraceJacToMap", 10);
966 }
967 
968 template <typename DataType, typename TypeNekBlkMatSharedPtr>
971  TypeNekBlkMatSharedPtr &gmtVar, const DataType &tmpDataType)
972 {
973  boost::ignore_unused(tmpDataType);
974 
975  int n1d = gmtxarray.size();
976  int n2d = gmtxarray[0].size();
977  int nConvectiveFields = n1d;
978 
979  ASSERTL0(n1d == n2d, "ElmtVarInvMtrx requires n1d==n2d");
980 
981  Array<OneD, unsigned int> rowSizes;
982  Array<OneD, unsigned int> colSizes;
983 
984  gmtxarray[0][0]->GetBlockSizes(rowSizes, colSizes);
985  int nTotElmt = rowSizes.size();
986  int nElmtCoef = rowSizes[0] - 1;
987  int nElmtCoef0 = -1;
988  int blocksize = -1;
989 
990  Array<OneD, unsigned int> tmprow(1);
991  TypeNekBlkMatSharedPtr tmpGmtx;
992 
993  Array<OneD, DataType> GMatData, ElmtMatData;
994  Array<OneD, DataType> tmpArray1, tmpArray2;
995 
996  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
997  {
998  int nrows = gmtxarray[0][0]->GetBlock(nelmt, nelmt)->GetRows();
999  int ncols = gmtxarray[0][0]->GetBlock(nelmt, nelmt)->GetColumns();
1000  ASSERTL0(nrows == ncols, "ElmtVarInvMtrx requires nrows==ncols");
1001 
1002  nElmtCoef = nrows;
1003 
1004  if (nElmtCoef0 != nElmtCoef)
1005  {
1006  nElmtCoef0 = nElmtCoef;
1007  int nElmtCoefVar = nElmtCoef0 * nConvectiveFields;
1008  blocksize = nElmtCoefVar * nElmtCoefVar;
1009  tmprow[0] = nElmtCoefVar;
1010  AllocateNekBlkMatDig(tmpGmtx, tmprow, tmprow);
1011  GMatData = tmpGmtx->GetBlock(0, 0)->GetPtr();
1012  }
1013 
1014  for (int n = 0; n < nConvectiveFields; n++)
1015  {
1016  for (int m = 0; m < nConvectiveFields; m++)
1017  {
1018  ElmtMatData = gmtxarray[m][n]->GetBlock(nelmt, nelmt)->GetPtr();
1019 
1020  for (int ncl = 0; ncl < nElmtCoef; ncl++)
1021  {
1022  int Goffset =
1023  (n * nElmtCoef + ncl) * nConvectiveFields * nElmtCoef +
1024  m * nElmtCoef;
1025  int Eoffset = ncl * nElmtCoef;
1026 
1027  Vmath::Vcopy(nElmtCoef, tmpArray1 = ElmtMatData + Eoffset,
1028  1, tmpArray2 = GMatData + Goffset, 1);
1029  }
1030  }
1031  }
1032 
1033  tmpGmtx->GetBlock(0, 0)->Invert();
1034 
1035  for (int m = 0; m < nConvectiveFields; m++)
1036  {
1037  for (int n = 0; n < nConvectiveFields; n++)
1038  {
1039  ElmtMatData = gmtxarray[m][n]->GetBlock(nelmt, nelmt)->GetPtr();
1040 
1041  for (int ncl = 0; ncl < nElmtCoef; ncl++)
1042  {
1043  int Goffset =
1044  (n * nElmtCoef + ncl) * nConvectiveFields * nElmtCoef +
1045  m * nElmtCoef;
1046  int Eoffset = ncl * nElmtCoef;
1047 
1048  Vmath::Vcopy(nElmtCoef, tmpArray1 = GMatData + Goffset, 1,
1049  tmpArray2 = ElmtMatData + Eoffset, 1);
1050  }
1051  }
1052  }
1053  ElmtMatData = gmtVar->GetBlock(nelmt, nelmt)->GetPtr();
1054  Vmath::Vcopy(blocksize, GMatData, 1, ElmtMatData, 1);
1055  }
1056  return;
1057 }
1058 
1059 template <typename DataType, typename TypeNekBlkMatSharedPtr>
1061  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1064  Array<OneD, TypeNekBlkMatSharedPtr> &TraceJacDeriv,
1065  Array<OneD, Array<OneD, DataType>> &TraceJacDerivSign,
1066  TensorOfArray5D<DataType> &TraceIPSymJacArray)
1067 {
1068  boost::ignore_unused(TraceJacDeriv, TraceJacDerivSign);
1069 
1070  int nvariables = inarray.size();
1071  int nTracePts = GetTraceTotPoints();
1072 
1073  // Store forwards/backwards space along trace space
1074  Array<OneD, Array<OneD, NekDouble>> Fwd(nvariables);
1075  Array<OneD, Array<OneD, NekDouble>> Bwd(nvariables);
1076 
1077  TypeNekBlkMatSharedPtr FJac, BJac;
1078  Array<OneD, unsigned int> n_blks1(nTracePts, nvariables);
1079 
1080  if (TraceJac.size() > 0)
1081  {
1082  FJac = TraceJac[0];
1083  BJac = TraceJac[1];
1084  }
1085  else
1086  {
1088 
1089  AllocateNekBlkMatDig(FJac, n_blks1, n_blks1);
1090  AllocateNekBlkMatDig(BJac, n_blks1, n_blks1);
1091  }
1092 
1094  {
1097  }
1098  else
1099  {
1100  for (int i = 0; i < nvariables; ++i)
1101  {
1102  Fwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1103  Bwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1104  m_fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
1105  }
1106  }
1107 
1109 
1110  NumCalcRiemFluxJac(nvariables, m_fields, AdvVel, inarray, qfield,
1111  m_bndEvaluateTime, Fwd, Bwd, FJac, BJac,
1112  TraceIPSymJacArray);
1113 
1114  TraceJac[0] = FJac;
1115  TraceJac[1] = BJac;
1116 }
1117 
1118 template <typename DataType, typename TypeNekBlkMatSharedPtr>
1120  const int nConvectiveFields,
1122  const Array<OneD, const Array<OneD, NekDouble>> &AdvVel,
1123  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1124  TensorOfArray3D<NekDouble> &qfield, const NekDouble &time,
1125  const Array<OneD, const Array<OneD, NekDouble>> &Fwd,
1126  const Array<OneD, const Array<OneD, NekDouble>> &Bwd,
1127  TypeNekBlkMatSharedPtr &FJac, TypeNekBlkMatSharedPtr &BJac,
1128  TensorOfArray5D<DataType> &TraceIPSymJacArray)
1129 {
1130  boost::ignore_unused(TraceIPSymJacArray);
1131 
1132  const NekDouble PenaltyFactor2 = 0.0;
1133  int nvariables = nConvectiveFields;
1134  int npoints = GetNpoints();
1135  // int nPts = npoints;
1136  int nTracePts = GetTraceTotPoints();
1137  int nDim = m_spacedim;
1138 
1139  Array<OneD, int> nonZeroIndex;
1140 
1141  Array<OneD, Array<OneD, NekDouble>> tmpinarry(nvariables);
1142  for (int i = 0; i < nvariables; i++)
1143  {
1144  tmpinarry[i] = Array<OneD, NekDouble>(npoints, 0.0);
1145  Vmath::Vcopy(npoints, inarray[i], 1, tmpinarry[i], 1);
1146  }
1147 
1148  // DmuDT of artificial diffusion is neglected
1149  // TODO: to consider the Jacobian of AV seperately
1152 
1153  Array<OneD, Array<OneD, NekDouble>> numflux(nvariables);
1154  for (int i = 0; i < nvariables; ++i)
1155  {
1156  numflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1157  }
1158 
1159  const MultiRegions::AssemblyMapDGSharedPtr TraceMap =
1160  fields[0]->GetTraceMap();
1161  TensorOfArray3D<NekDouble> qBwd(nDim);
1162  TensorOfArray3D<NekDouble> qFwd(nDim);
1163  if (m_viscousJacFlag)
1164  {
1165  for (int nd = 0; nd < nDim; ++nd)
1166  {
1167  qBwd[nd] = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
1168  qFwd[nd] = Array<OneD, Array<OneD, NekDouble>>(nConvectiveFields);
1169  for (int i = 0; i < nConvectiveFields; ++i)
1170  {
1171  qBwd[nd][i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1172  qFwd[nd][i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1173 
1174  fields[i]->GetFwdBwdTracePhys(qfield[nd][i], qFwd[nd][i],
1175  qBwd[nd][i], true, true, false);
1176  TraceMap->GetAssemblyCommDG()->PerformExchange(qFwd[nd][i],
1177  qBwd[nd][i]);
1178  }
1179  }
1180  }
1181 
1182  CalcTraceNumericalFlux(nConvectiveFields, nDim, npoints, nTracePts,
1183  PenaltyFactor2, fields, AdvVel, inarray, time,
1184  qfield, Fwd, Bwd, qFwd, qBwd, MuVarTrace,
1185  nonZeroIndex, numflux);
1186 
1187  int nFields = nvariables;
1188  Array<OneD, Array<OneD, NekDouble>> plusFwd(nFields), plusBwd(nFields);
1189  Array<OneD, Array<OneD, NekDouble>> Jacvect(nFields);
1190  Array<OneD, Array<OneD, NekDouble>> FwdBnd(nFields);
1191  Array<OneD, Array<OneD, NekDouble>> plusflux(nFields);
1192  for (int i = 0; i < nFields; i++)
1193  {
1194  Jacvect[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1195  plusFwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1196  plusBwd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1197  plusflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1198  FwdBnd[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1199  }
1200 
1201  for (int i = 0; i < nFields; i++)
1202  {
1203  Vmath::Vcopy(nTracePts, Fwd[i], 1, plusFwd[i], 1);
1204  Vmath::Vcopy(nTracePts, Bwd[i], 1, plusBwd[i], 1);
1205  }
1206 
1207  NekDouble eps = 1.0E-6;
1208 
1209  Array<OneD, DataType> tmpMatData;
1210  // Fwd Jacobian
1211  for (int i = 0; i < nFields; i++)
1212  {
1213  NekDouble epsvar = eps * m_magnitdEstimat[i];
1214  NekDouble oepsvar = 1.0 / epsvar;
1215  Vmath::Sadd(nTracePts, epsvar, Fwd[i], 1, plusFwd[i], 1);
1216 
1217  if (m_bndConds.size())
1218  {
1219  for (int i = 0; i < nFields; i++)
1220  {
1221  Vmath::Vcopy(nTracePts, plusFwd[i], 1, FwdBnd[i], 1);
1222  }
1223  // Loop over user-defined boundary conditions
1224  for (auto &x : m_bndConds)
1225  {
1226  x->Apply(FwdBnd, tmpinarry, time);
1227  }
1228  }
1229 
1230  for (int j = 0; j < nFields; j++)
1231  {
1232  m_fields[j]->FillBwdWithBoundCond(plusFwd[j], plusBwd[j]);
1233  }
1234 
1235  CalcTraceNumericalFlux(nConvectiveFields, nDim, npoints, nTracePts,
1236  PenaltyFactor2, fields, AdvVel, inarray, time,
1237  qfield, plusFwd, plusBwd, qFwd, qBwd, MuVarTrace,
1238  nonZeroIndex, plusflux);
1239 
1240  for (int n = 0; n < nFields; n++)
1241  {
1242  Vmath::Vsub(nTracePts, plusflux[n], 1, numflux[n], 1, Jacvect[n],
1243  1);
1244  Vmath::Smul(nTracePts, oepsvar, Jacvect[n], 1, Jacvect[n], 1);
1245  }
1246  for (int j = 0; j < nTracePts; j++)
1247  {
1248  tmpMatData = FJac->GetBlock(j, j)->GetPtr();
1249  for (int n = 0; n < nFields; n++)
1250  {
1251  tmpMatData[n + i * nFields] = DataType(Jacvect[n][j]);
1252  }
1253  }
1254 
1255  Vmath::Vcopy(nTracePts, Fwd[i], 1, plusFwd[i], 1);
1256  }
1257 
1258  // Reset the boundary conditions
1259  if (m_bndConds.size())
1260  {
1261  for (int i = 0; i < nFields; i++)
1262  {
1263  Vmath::Vcopy(nTracePts, Fwd[i], 1, FwdBnd[i], 1);
1264  }
1265  // Loop over user-defined boundary conditions
1266  for (auto &x : m_bndConds)
1267  {
1268  x->Apply(FwdBnd, tmpinarry, time);
1269  }
1270  }
1271 
1272  for (int i = 0; i < nFields; i++)
1273  {
1274  Vmath::Vcopy(nTracePts, Bwd[i], 1, plusBwd[i], 1);
1275  }
1276 
1277  for (int i = 0; i < nFields; i++)
1278  {
1279  NekDouble epsvar = eps * m_magnitdEstimat[i];
1280  NekDouble oepsvar = 1.0 / epsvar;
1281 
1282  Vmath::Sadd(nTracePts, epsvar, Bwd[i], 1, plusBwd[i], 1);
1283 
1284  for (int j = 0; j < nFields; j++)
1285  {
1286  m_fields[j]->FillBwdWithBoundCond(Fwd[j], plusBwd[j]);
1287  }
1288 
1289  CalcTraceNumericalFlux(nConvectiveFields, nDim, npoints, nTracePts,
1290  PenaltyFactor2, fields, AdvVel, inarray, time,
1291  qfield, Fwd, plusBwd, qFwd, qBwd, MuVarTrace,
1292  nonZeroIndex, plusflux);
1293 
1294  for (int n = 0; n < nFields; n++)
1295  {
1296  Vmath::Vsub(nTracePts, plusflux[n], 1, numflux[n], 1, Jacvect[n],
1297  1);
1298  Vmath::Smul(nTracePts, oepsvar, Jacvect[n], 1, Jacvect[n], 1);
1299  }
1300  for (int j = 0; j < nTracePts; j++)
1301  {
1302  tmpMatData = BJac->GetBlock(j, j)->GetPtr();
1303  for (int n = 0; n < nFields; n++)
1304  {
1305  tmpMatData[n + i * nFields] = DataType(Jacvect[n][j]);
1306  }
1307  }
1308 
1309  Vmath::Vcopy(nTracePts, Bwd[i], 1, plusBwd[i], 1);
1310  }
1311 }
1312 
1314  const int nConvectiveFields, const int nDim, const int nPts,
1315  const int nTracePts, const NekDouble PenaltyFactor2,
1317  const Array<OneD, const Array<OneD, NekDouble>> &AdvVel,
1318  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1319  const NekDouble time, TensorOfArray3D<NekDouble> &qfield,
1320  const Array<OneD, const Array<OneD, NekDouble>> &vFwd,
1321  const Array<OneD, const Array<OneD, NekDouble>> &vBwd,
1322  const Array<OneD, const TensorOfArray2D<NekDouble>> &qFwd,
1323  const Array<OneD, const TensorOfArray2D<NekDouble>> &qBwd,
1324  const Array<OneD, NekDouble> &MuVarTrace, Array<OneD, int> &nonZeroIndex,
1325  Array<OneD, Array<OneD, NekDouble>> &traceflux)
1326 {
1327  boost::ignore_unused(nDim, nPts, PenaltyFactor2, time, qFwd, qBwd,
1328  MuVarTrace);
1329 
1330  if (m_advectionJacFlag)
1331  {
1332  m_advObject->AdvectTraceFlux(nConvectiveFields, m_fields, AdvVel,
1333  inarray, traceflux, m_bndEvaluateTime,
1334  vFwd, vBwd);
1335  }
1336  else
1337  {
1338  for (int i = 0; i < nConvectiveFields; i++)
1339  {
1340  traceflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1341  }
1342  }
1343 
1344  if (m_viscousJacFlag)
1345  {
1346  Array<OneD, Array<OneD, NekDouble>> visflux(nConvectiveFields);
1347  for (int i = 0; i < nConvectiveFields; i++)
1348  {
1349  visflux[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
1350  }
1351 
1352  m_diffusion->DiffuseTraceFlux(fields, inarray, qfield,
1354  vFwd, vBwd, nonZeroIndex);
1355  for (int i = 0; i < nConvectiveFields; i++)
1356  {
1357  Vmath::Vsub(nTracePts, traceflux[i], 1, visflux[i], 1, traceflux[i],
1358  1);
1359  }
1360  }
1361 }
1362 
1364  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1367  Array<OneD, SNekBlkMatSharedPtr> &TraceJacDeriv,
1368  Array<OneD, Array<OneD, NekSingle>> &TraceJacDerivSign,
1369  TensorOfArray4D<NekSingle> &TraceJacArray,
1370  TensorOfArray4D<NekSingle> &TraceJacDerivArray,
1371  TensorOfArray5D<NekSingle> &TraceIPSymJacArray)
1372 {
1374 
1375  if (m_viscousJacFlag)
1376  {
1377  CalcPhysDeriv(inarray, qfield);
1378  }
1379 
1380  NekSingle zero = 0.0;
1381  Fill2DArrayOfBlkDiagonalMat(gmtxarray, zero);
1382 
1383  LibUtilities::Timer timer;
1384  timer.Start();
1385  AddMatNSBlkDiagVol(inarray, qfield, gmtxarray, m_stdSMatDataDBB,
1387  timer.Stop();
1388  timer.AccumulateRegion("CFSImplicit::AddMatNSBlkDiagVol", 2);
1389 
1390  timer.Start();
1391  AddMatNSBlkDiagBnd(inarray, qfield, gmtxarray, TraceJac, TraceJacDeriv,
1392  TraceJacDerivSign, TraceIPSymJacArray);
1393  timer.Stop();
1394  timer.AccumulateRegion("CFSImplicit::AddMatNSBlkDiagBnd", 2);
1395 
1397 
1398  timer.Start();
1399  ElmtVarInvMtrx(gmtxarray, gmtVar, zero);
1400  timer.Stop();
1401  timer.AccumulateRegion("CFSImplicit::ElmtVarInvMtrx", 2);
1402 
1403  TransTraceJacMatToArray(TraceJac, TraceJacDeriv, TraceJacArray,
1404  TraceJacDerivArray);
1405 }
1406 
1408  const int nConvectiveFields, const int nElmtPnt,
1409  const Array<OneD, const Array<OneD, NekDouble>> &locVars,
1410  const TensorOfArray3D<NekDouble> &locDerv,
1411  const Array<OneD, NekDouble> &locmu, const Array<OneD, NekDouble> &locDmuDT,
1412  const Array<OneD, NekDouble> &normals, DNekMatSharedPtr &wspMat,
1413  Array<OneD, Array<OneD, NekDouble>> &PntJacArray)
1414 {
1415  boost::ignore_unused(nConvectiveFields, nElmtPnt, locVars, locDerv, locmu,
1416  locDmuDT, normals, wspMat, PntJacArray);
1417  // Do nothing by default
1418 }
1419 
1420 template <typename DataType, typename TypeNekBlkMatSharedPtr>
1423  const NekDouble dtlamda, const DataType tmpDataType)
1424 {
1425  boost::ignore_unused(tmpDataType);
1426 
1428  std::shared_ptr<LocalRegions::ExpansionVector> pexp = explist->GetExp();
1429  int nTotElmt = (*pexp).size();
1430  int nConvectiveFields = m_fields.size();
1431 
1432  NekDouble Negdtlamda = -dtlamda;
1433 
1434  Array<OneD, NekDouble> pseudotimefactor(nTotElmt, 0.0);
1435  Vmath::Fill(nTotElmt, Negdtlamda, pseudotimefactor, 1);
1436 
1437  Array<OneD, DataType> GMatData;
1438  for (int m = 0; m < nConvectiveFields; m++)
1439  {
1440  for (int n = 0; n < nConvectiveFields; n++)
1441  {
1442  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
1443  {
1444  GMatData = gmtxarray[m][n]->GetBlock(nelmt, nelmt)->GetPtr();
1445  DataType factor = DataType(pseudotimefactor[nelmt]);
1446 
1447  Vmath::Smul(GMatData.size(), factor, GMatData, 1, GMatData, 1);
1448  }
1449  }
1450  }
1451 
1452  DNekMatSharedPtr MassMat;
1453  Array<OneD, NekDouble> BwdMatData, MassMatData, tmp;
1455  Array<OneD, DataType> MassMatDataDataType;
1457 
1458  for (int nelmt = 0; nelmt < nTotElmt; nelmt++)
1459  {
1460  int nelmtcoef = GetNcoeffs(nelmt);
1461  int nelmtpnts = GetTotPoints(nelmt);
1462  LibUtilities::ShapeType ElmtTypeNow =
1463  explist->GetExp(nelmt)->DetShapeType();
1464 
1465  if (tmp.size() != nelmtcoef || (ElmtTypeNow != ElmtTypePrevious))
1466  {
1468  stdExp = explist->GetExp(nelmt)->GetStdExp();
1470  stdExp->DetShapeType(), *stdExp);
1471 
1472  DNekMatSharedPtr BwdMat = stdExp->GetStdMatrix(matkey);
1473  BwdMatData = BwdMat->GetPtr();
1474 
1475  if (nelmtcoef != tmp.size())
1476  {
1477  tmp = Array<OneD, NekDouble>(nelmtcoef, 0.0);
1479  nelmtcoef, nelmtcoef, 0.0);
1480  MassMatData = MassMat->GetPtr();
1481  MassMatDataDataType =
1482  Array<OneD, DataType>(nelmtcoef * nelmtcoef);
1483  }
1484 
1485  ElmtTypePrevious = ElmtTypeNow;
1486  }
1487 
1488  for (int np = 0; np < nelmtcoef; np++)
1489  {
1490  explist->GetExp(nelmt)->IProductWRTBase(BwdMatData + np * nelmtpnts,
1491  tmp);
1492  Vmath::Vcopy(nelmtcoef, tmp, 1, tmp2 = MassMatData + np * nelmtcoef,
1493  1);
1494  }
1495  for (int i = 0; i < MassMatData.size(); i++)
1496  {
1497  MassMatDataDataType[i] = DataType(MassMatData[i]);
1498  }
1499 
1500  for (int m = 0; m < nConvectiveFields; m++)
1501  {
1502  GMatData = gmtxarray[m][m]->GetBlock(nelmt, nelmt)->GetPtr();
1503  Vmath::Vadd(MassMatData.size(), MassMatDataDataType, 1, GMatData, 1,
1504  GMatData, 1);
1505  }
1506  }
1507  return;
1508 }
1509 
1510 template <typename DataType, typename TypeNekBlkMatSharedPtr>
1512  const Array<OneD, TypeNekBlkMatSharedPtr> &TraceJac,
1513  const Array<OneD, TypeNekBlkMatSharedPtr> &TraceJacDeriv,
1514  TensorOfArray4D<DataType> &TraceJacArray,
1515  TensorOfArray4D<DataType> &TraceJacDerivArray)
1516 {
1517  boost::ignore_unused(TraceJacArray, TraceJacDeriv, TraceJacDerivArray);
1518 
1519  int nFwdBwd, nDiagBlks, nvar0Jac, nvar1Jac;
1520 
1521  Array<OneD, unsigned int> rowSizes;
1522  Array<OneD, unsigned int> colSizes;
1523  nFwdBwd = TraceJac.size();
1524  TraceJac[0]->GetBlockSizes(rowSizes, colSizes);
1525  nDiagBlks = rowSizes.size();
1526  nvar0Jac = rowSizes[1] - rowSizes[0];
1527  nvar1Jac = colSizes[1] - colSizes[0];
1528 
1529  if (0 == TraceJacArray.size())
1530  {
1531  TraceJacArray = TensorOfArray4D<DataType>(nFwdBwd);
1532  for (int nlr = 0; nlr < nFwdBwd; nlr++)
1533  {
1534  TraceJacArray[nlr] = TensorOfArray3D<DataType>(nvar0Jac);
1535  for (int m = 0; m < nvar0Jac; m++)
1536  {
1537  TraceJacArray[nlr][m] =
1539  for (int n = 0; n < nvar1Jac; n++)
1540  {
1541  TraceJacArray[nlr][m][n] = Array<OneD, DataType>(nDiagBlks);
1542  }
1543  }
1544  }
1545  }
1546 
1547  for (int nlr = 0; nlr < nFwdBwd; nlr++)
1548  {
1549  const TypeNekBlkMatSharedPtr tmpMat = TraceJac[nlr];
1550  TensorOfArray3D<DataType> tmpaa = TraceJacArray[nlr];
1551  TranSamesizeBlkDiagMatIntoArray(tmpMat, tmpaa);
1552  }
1553 
1554  return;
1555 }
1556 
1558  const MultiRegions::ExpListSharedPtr &explist,
1559  const Array<OneD, const Array<OneD, NekDouble>> &normals,
1560  const int nDervDir,
1561  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1562  TensorOfArray5D<NekDouble> &ElmtJacArray, const int nFluxDir)
1563 {
1564  boost::ignore_unused(explist, normals, nDervDir, inarray, ElmtJacArray,
1565  nFluxDir);
1566  NEKERROR(ErrorUtil::efatal, "v_GetFluxDerivJacDirctn not coded");
1567 }
1568 
1570  const int nConvectiveFields, const int nElmtPnt, const int nDervDir,
1571  const Array<OneD, const Array<OneD, NekDouble>> &locVars,
1572  const Array<OneD, NekDouble> &locmu,
1573  const Array<OneD, const Array<OneD, NekDouble>> &locnormal,
1574  DNekMatSharedPtr &wspMat, Array<OneD, Array<OneD, NekDouble>> &PntJacArray)
1575 {
1576  boost::ignore_unused(nConvectiveFields, nElmtPnt, nDervDir, locVars, locmu,
1577  locnormal, wspMat, PntJacArray);
1578  NEKERROR(ErrorUtil::efatal, "v_GetFluxDerivJacDirctn not coded");
1579 }
1580 
1582  const MultiRegions::ExpListSharedPtr &explist,
1583  const Array<OneD, const Array<OneD, NekDouble>> &normals,
1584  const int nDervDir,
1585  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1587 {
1588  boost::ignore_unused(explist, normals, nDervDir, inarray, ElmtJac);
1589 }
1590 
1592  const int nConvectiveFields, const int nElmtPnt,
1593  const Array<OneD, const Array<OneD, NekDouble>> &locVars,
1594  const Array<OneD, NekDouble> &normals, DNekMatSharedPtr &wspMat,
1595  Array<OneD, Array<OneD, NekDouble>> &PntJacArray)
1596 {
1597  Array<OneD, NekDouble> wspMatData = wspMat->GetPtr();
1598 
1599  int matsize = nConvectiveFields * nConvectiveFields;
1600 
1601  Array<OneD, NekDouble> pointVar(nConvectiveFields);
1602 
1603  for (int npnt = 0; npnt < nElmtPnt; npnt++)
1604  {
1605  for (int j = 0; j < nConvectiveFields; j++)
1606  {
1607  pointVar[j] = locVars[j][npnt];
1608  }
1609 
1610  GetFluxVectorJacPoint(nConvectiveFields, pointVar, normals, wspMat);
1611 
1612  Vmath::Vcopy(matsize, wspMatData, 1, PntJacArray[npnt], 1);
1613  }
1614  return;
1615 }
1616 
1618  const int nConvectiveFields, const Array<OneD, NekDouble> &conservVar,
1619  const Array<OneD, NekDouble> &normals, DNekMatSharedPtr &fluxJac)
1620 {
1621  int nvariables = conservVar.size();
1622  const int nvariables3D = 5;
1623  int expDim = m_spacedim;
1624 
1625  NekDouble fsw, efix_StegerWarming;
1626  efix_StegerWarming = 0.0;
1627  fsw = 0.0; // exact flux Jacobian if fsw=0.0
1628  if (nvariables > expDim + 2)
1629  {
1630  NEKERROR(ErrorUtil::efatal, "nvariables > expDim+2 case not coded")
1631  }
1632 
1633  Array<OneD, NekDouble> fluxJacData;
1634  ;
1635  fluxJacData = fluxJac->GetPtr();
1636 
1637  if (nConvectiveFields == nvariables3D)
1638  {
1639  PointFluxJacobianPoint(conservVar, normals, fluxJac, efix_StegerWarming,
1640  fsw);
1641  }
1642  else
1643  {
1644  DNekMatSharedPtr PointFJac3D =
1646  nvariables3D, 0.0);
1647 
1648  Array<OneD, NekDouble> PointFJac3DData;
1649  PointFJac3DData = PointFJac3D->GetPtr();
1650 
1651  Array<OneD, NekDouble> PointFwd(nvariables3D, 0.0);
1652 
1653  Array<OneD, unsigned int> index(nvariables);
1654 
1655  index[nvariables - 1] = 4;
1656  for (int i = 0; i < nvariables - 1; i++)
1657  {
1658  index[i] = i;
1659  }
1660 
1661  int nj = 0;
1662  int nk = 0;
1663  for (int j = 0; j < nvariables; j++)
1664  {
1665  nj = index[j];
1666  PointFwd[nj] = conservVar[j];
1667  }
1668 
1669  PointFluxJacobianPoint(PointFwd, normals, PointFJac3D,
1670  efix_StegerWarming, fsw);
1671 
1672  for (int j = 0; j < nvariables; j++)
1673  {
1674  nj = index[j];
1675  for (int k = 0; k < nvariables; k++)
1676  {
1677  nk = index[k];
1678  fluxJacData[j + k * nConvectiveFields] =
1679  PointFJac3DData[nj + nk * nvariables3D];
1680  }
1681  }
1682  }
1683 }
1684 
1685 /**
1686  * @brief Add the diffusions terms to the right-hand side
1687  * Similar to DoDiffusion() but with outarray in coefficient space
1688  */
1690  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
1691  Array<OneD, Array<OneD, NekDouble>> &outarray,
1692  const Array<OneD, const Array<OneD, NekDouble>> &pFwd,
1693  const Array<OneD, const Array<OneD, NekDouble>> &pBwd)
1694 {
1695  v_DoDiffusionCoeff(inarray, outarray, pFwd, pBwd);
1696 }
1697 
1700  const bool &flag)
1701 {
1702  boost::ignore_unused(flag);
1703  const Array<OneD, const NekDouble> refsol = m_nonlinsol->GetRefSolution();
1704  const Array<OneD, const NekDouble> refres = m_nonlinsol->GetRefResidual();
1705  const Array<OneD, const NekDouble> refsource =
1706  m_nonlinsol->GetRefSourceVec();
1707 
1708  NekDouble eps = m_jacobiFreeEps;
1709  unsigned int nTotalGlobal = inarray.size();
1710  NekDouble magninarray = Vmath::Dot(nTotalGlobal, inarray, inarray);
1711  m_comm->AllReduce(magninarray, Nektar::LibUtilities::ReduceSum);
1712  eps *= sqrt((sqrt(m_inArrayNorm) + 1.0) / magninarray);
1713 
1714  NekDouble oeps = 1.0 / eps;
1715  unsigned int ntotal = inarray.size();
1716  Array<OneD, NekDouble> solplus{ntotal, 0.0};
1717  Array<OneD, NekDouble> resplus{ntotal, 0.0};
1718 
1719  Vmath::Svtvp(ntotal, eps, inarray, 1, refsol, 1, solplus, 1);
1720 
1721  NonlinSysEvaluatorCoeff(solplus, resplus, flag, refsource);
1722 
1723  Vmath::Vsub(ntotal, resplus, 1, refres, 1, out, 1);
1724  Vmath::Smul(ntotal, oeps, out, 1, out, 1);
1725 
1726  return;
1727 }
1728 
1729 template <typename DataType, typename TypeNekBlkMatSharedPtr>
1731  const TypeNekBlkMatSharedPtr &BlkMat, TensorOfArray3D<DataType> &MatArray)
1732 {
1733  Array<OneD, unsigned int> rowSizes;
1734  Array<OneD, unsigned int> colSizes;
1735  BlkMat->GetBlockSizes(rowSizes, colSizes);
1736  int nDiagBlks = rowSizes.size();
1737  int nvar0 = rowSizes[1] - rowSizes[0];
1738  int nvar1 = colSizes[1] - colSizes[0];
1739 
1740  Array<OneD, DataType> ElmtMatData;
1741 
1742  for (int i = 0; i < nDiagBlks; i++)
1743  {
1744  ElmtMatData = BlkMat->GetBlock(i, i)->GetPtr();
1745  for (int n = 0; n < nvar1; n++)
1746  {
1747  int noffset = n * nvar0;
1748  for (int m = 0; m < nvar0; m++)
1749  {
1750  MatArray[m][n][i] = ElmtMatData[m + noffset];
1751  }
1752  }
1753  }
1754 }
1755 
1756 template <typename DataType, typename TypeNekBlkMatSharedPtr>
1759  const DataType valu)
1760 {
1761 
1762  int n1d = gmtxarray.size();
1763 
1764  for (int n1 = 0; n1 < n1d; ++n1)
1765  {
1766  Fill1DArrayOfBlkDiagonalMat(gmtxarray[n1], valu);
1767  }
1768 }
1769 
1770 template <typename DataType, typename TypeNekBlkMatSharedPtr>
1772  Array<OneD, TypeNekBlkMatSharedPtr> &gmtxarray, const DataType valu)
1773 {
1774  int n1d = gmtxarray.size();
1775 
1776  Array<OneD, unsigned int> rowSizes;
1777  Array<OneD, unsigned int> colSizes;
1778 
1779  Array<OneD, DataType> loc_mat_arr;
1780 
1781  for (int n1 = 0; n1 < n1d; ++n1)
1782  {
1783  gmtxarray[n1]->GetBlockSizes(rowSizes, colSizes);
1784  int nelmts = rowSizes.size();
1785 
1786  for (int i = 0; i < nelmts; ++i)
1787  {
1788  loc_mat_arr = gmtxarray[n1]->GetBlock(i, i)->GetPtr();
1789 
1790  int nrows = gmtxarray[n1]->GetBlock(i, i)->GetRows();
1791  int ncols = gmtxarray[n1]->GetBlock(i, i)->GetColumns();
1792 
1793  Vmath::Fill(nrows * ncols, valu, loc_mat_arr, 1);
1794  }
1795  }
1796 }
1797 
1798 // Currently duplacate in compressibleFlowSys
1799 // if fsw=+-1 calculate the steger-Warming flux vector splitting flux Jacobian
1800 // if fsw=0 calculate the Jacobian of the exact flux
1801 // efix is the numerical flux entropy fix parameter
1803  const Array<OneD, NekDouble> &normals,
1804  DNekMatSharedPtr &FJac,
1805  const NekDouble efix,
1806  const NekDouble fsw)
1807 {
1808  Array<OneD, NekDouble> FJacData = FJac->GetPtr();
1809  const int nvariables3D = 5;
1810 
1811  NekDouble ro, vx, vy, vz, ps, gama, ae;
1812  NekDouble a, a2, h, h0, v2, vn, eps, eps2;
1813  NekDouble nx, ny, nz;
1814  NekDouble sn, osn, nxa, nya, nza, vna;
1815  NekDouble l1, l4, l5, al1, al4, al5, x1, x2, x3, y1;
1816  NekDouble c1, d1, c2, d2, c3, d3, c4, d4, c5, d5;
1817  NekDouble sml_ssf = 1.0E-12;
1818 
1819  NekDouble fExactorSplt = 2.0 - abs(fsw); // if fsw=+-1 calculate
1820 
1821  NekDouble rhoL = Fwd[0];
1822  NekDouble rhouL = Fwd[1];
1823  NekDouble rhovL = Fwd[2];
1824  NekDouble rhowL = Fwd[3];
1825  NekDouble EL = Fwd[4];
1826 
1827  ro = rhoL;
1828  vx = rhouL / rhoL;
1829  vy = rhovL / rhoL;
1830  vz = rhowL / rhoL;
1831 
1832  // Internal energy (per unit mass)
1833  NekDouble eL = (EL - 0.5 * (rhouL * vx + rhovL * vy + rhowL * vz)) / rhoL;
1834 
1835  ps = m_varConv->Geteos()->GetPressure(rhoL, eL);
1836  gama = m_gamma;
1837 
1838  ae = gama - 1.0;
1839  v2 = vx * vx + vy * vy + vz * vz;
1840  a2 = gama * ps / ro;
1841  h = a2 / ae;
1842 
1843  h0 = h + 0.5 * v2;
1844  a = sqrt(a2);
1845 
1846  nx = normals[0];
1847  ny = normals[1];
1848  nz = normals[2];
1849  vn = nx * vx + ny * vy + nz * vz;
1850  sn = std::max(sqrt(nx * nx + ny * ny + nz * nz), sml_ssf);
1851  osn = 1.0 / sn;
1852 
1853  nxa = nx * osn;
1854  nya = ny * osn;
1855  nza = nz * osn;
1856  vna = vn * osn;
1857  l1 = vn;
1858  l4 = vn + sn * a;
1859  l5 = vn - sn * a;
1860 
1861  eps = efix * sn;
1862  eps2 = eps * eps;
1863 
1864  al1 = sqrt(l1 * l1 + eps2);
1865  al4 = sqrt(l4 * l4 + eps2);
1866  al5 = sqrt(l5 * l5 + eps2);
1867 
1868  l1 = 0.5 * (fExactorSplt * l1 + fsw * al1);
1869  l4 = 0.5 * (fExactorSplt * l4 + fsw * al4);
1870  l5 = 0.5 * (fExactorSplt * l5 + fsw * al5);
1871 
1872  x1 = 0.5 * (l4 + l5);
1873  x2 = 0.5 * (l4 - l5);
1874  x3 = x1 - l1;
1875  y1 = 0.5 * v2;
1876  c1 = ae * x3 / a2;
1877  d1 = x2 / a;
1878 
1879  int nVar0 = 0;
1880  int nVar1 = nvariables3D;
1881  int nVar2 = 2 * nvariables3D;
1882  int nVar3 = 3 * nvariables3D;
1883  int nVar4 = 4 * nvariables3D;
1884  FJacData[nVar0] = c1 * y1 - d1 * vna + l1;
1885  FJacData[nVar1] = -c1 * vx + d1 * nxa;
1886  FJacData[nVar2] = -c1 * vy + d1 * nya;
1887  FJacData[nVar3] = -c1 * vz + d1 * nza;
1888  FJacData[nVar4] = c1;
1889  c2 = c1 * vx + d1 * nxa * ae;
1890  d2 = x3 * nxa + d1 * vx;
1891  FJacData[1 + nVar0] = c2 * y1 - d2 * vna;
1892  FJacData[1 + nVar1] = -c2 * vx + d2 * nxa + l1;
1893  FJacData[1 + nVar2] = -c2 * vy + d2 * nya;
1894  FJacData[1 + nVar3] = -c2 * vz + d2 * nza;
1895  FJacData[1 + nVar4] = c2;
1896  c3 = c1 * vy + d1 * nya * ae;
1897  d3 = x3 * nya + d1 * vy;
1898  FJacData[2 + nVar0] = c3 * y1 - d3 * vna;
1899  FJacData[2 + nVar1] = -c3 * vx + d3 * nxa;
1900  FJacData[2 + nVar2] = -c3 * vy + d3 * nya + l1;
1901  FJacData[2 + nVar3] = -c3 * vz + d3 * nza;
1902  FJacData[2 + nVar4] = c3;
1903  c4 = c1 * vz + d1 * nza * ae;
1904  d4 = x3 * nza + d1 * vz;
1905  FJacData[3 + nVar0] = c4 * y1 - d4 * vna;
1906  FJacData[3 + nVar1] = -c4 * vx + d4 * nxa;
1907  FJacData[3 + nVar2] = -c4 * vy + d4 * nya;
1908  FJacData[3 + nVar3] = -c4 * vz + d4 * nza + l1;
1909  FJacData[3 + nVar4] = c4;
1910  c5 = c1 * h0 + d1 * vna * ae;
1911  d5 = x3 * vna + d1 * h0;
1912  FJacData[4 + nVar0] = c5 * y1 - d5 * vna;
1913  FJacData[4 + nVar1] = -c5 * vx + d5 * nxa;
1914  FJacData[4 + nVar2] = -c5 * vy + d5 * nya;
1915  FJacData[4 + nVar3] = -c5 * vz + d5 * nza;
1916  FJacData[4 + nVar4] = c5 + l1;
1917 }
1918 
1920 {
1921  bool flag =
1922  (m_time + m_timestep > m_fintime && m_fintime > 0.0) ||
1924  if (m_explicitAdvection)
1925  {
1926  flag = true;
1927  }
1928  else
1929  {
1930  flag = flag || m_preconCfs->UpdatePreconMatCheck(NullNekDouble1DArray,
1932  }
1933  return flag;
1934 }
1935 
1936 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
size_type size() const
Returns the array's size.
void DoDiffusionCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, const Array< OneD, NekDouble >> &pFwd, const Array< OneD, const Array< OneD, NekDouble >> &pBwd)
Add the diffusions terms to the right-hand side Similar to DoDiffusion() but with outarray in coeffic...
void CalcVolJacStdMat(TensorOfArray4D< DataType > &StdMatDataDBB, TensorOfArray5D< DataType > &StdMatDataDBDB)
void ElmtVarInvMtrx(Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr >> &gmtxarray, TypeNekBlkMatSharedPtr &gmtVar, const DataType &tmpDatatype)
void DoAdvectionCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble >> &pFwd, const Array< OneD, const Array< OneD, NekDouble >> &pBwd)
Compute the advection terms for the right-hand side.
virtual void v_DoDiffusionCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const Array< OneD, const Array< OneD, NekDouble >> &pFwd, const Array< OneD, const Array< OneD, NekDouble >> &pBwd)
LibUtilities::NekNonlinSysSharedPtr m_nonlinsol
void AddMatNSBlkDiagVol(const Array< OneD, const Array< OneD, NekDouble >> &inarray, const Array< OneD, const TensorOfArray2D< NekDouble >> &qfield, Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr >> &gmtxarray, TensorOfArray4D< DataType > &StdMatDataDBB, TensorOfArray5D< DataType > &StdMatDataDBDB)
void PointFluxJacobianPoint(const Array< OneD, NekDouble > &Fwd, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &FJac, const NekDouble efix, const NekDouble fsw)
void PreconCoeff(const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const bool &flag)
void MultiplyElmtInvMassPlusSource(Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr >> &gmtxarray, const NekDouble dtlamda, const DataType tmpDataType)
void NonlinSysEvaluatorCoeff1D(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out, const bool &flag)
void CalcRefValues(const Array< OneD, const NekDouble > &inarray)
virtual void v_GetFluxDerivJacDirctn(const MultiRegions::ExpListSharedPtr &explist, const Array< OneD, const Array< OneD, NekDouble >> &normals, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray5D< NekDouble > &ElmtJacArray, const int nFluxDir)
void GetFluxVectorJacDirElmt(const int nConvectiveFields, const int nElmtPnt, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
void Fill2DArrayOfBlkDiagonalMat(Array< OneD, Array< OneD, TypeNekBlkMatSharedPtr >> &gmtxarray, const DataType valu)
void DoImplicitSolveCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inpnts, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out, const NekDouble time, const NekDouble lambda)
void CalcTraceNumericalFlux(const int nConvectiveFields, const int nDim, const int nPts, const int nTracePts, const NekDouble PenaltyFactor2, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const Array< OneD, NekDouble >> &AdvVel, const Array< OneD, const Array< OneD, NekDouble >> &inarray, const NekDouble time, TensorOfArray3D< NekDouble > &qfield, const Array< OneD, const Array< OneD, NekDouble >> &vFwd, const Array< OneD, const Array< OneD, NekDouble >> &vBwd, const Array< OneD, const TensorOfArray2D< NekDouble >> &qFwd, const Array< OneD, const TensorOfArray2D< NekDouble >> &qBwd, const Array< OneD, NekDouble > &MuVarTrace, Array< OneD, int > &nonZeroIndex, Array< OneD, Array< OneD, NekDouble >> &traceflux)
Array< OneD, Array< OneD, NekDouble > > m_solutionPhys
void AllocateNekBlkMatDig(SNekBlkMatSharedPtr &mat, const Array< OneD, unsigned int > nrow, const Array< OneD, unsigned int > ncol)
TensorOfArray4D< NekSingle > m_stdSMatDataDBB
void NonlinSysEvaluatorCoeff(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out, const bool &flag=true, const Array< OneD, const NekDouble > &source=NullNekDouble1DArray)
void TranSamesizeBlkDiagMatIntoArray(const TypeNekBlkMatSharedPtr &BlkMat, TensorOfArray3D< DataType > &MatArray)
void NumCalcRiemFluxJac(const int nConvectiveFields, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, const Array< OneD, NekDouble >> &AdvVel, const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield, const NekDouble &time, const Array< OneD, const Array< OneD, NekDouble >> &Fwd, const Array< OneD, const Array< OneD, NekDouble >> &Bwd, TypeNekBlkMatSharedPtr &FJac, TypeNekBlkMatSharedPtr &BJac, TensorOfArray5D< DataType > &TraceIPSymJacArray)
virtual ~CFSImplicit()
Destructor for CFSImplicit class.
void DoOdeRhsCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the right-hand side.
void GetTraceJac(const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield, Array< OneD, TypeNekBlkMatSharedPtr > &TraceJac, Array< OneD, TypeNekBlkMatSharedPtr > &TraceJacDeriv, Array< OneD, Array< OneD, DataType >> &TraceJacDerivSign, TensorOfArray5D< DataType > &TraceIPSymJacArray)
void GetFluxVectorJacPoint(const int nConvectiveFields, const Array< OneD, NekDouble > &conservVar, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &fluxJac)
TensorOfArray5D< NekSingle > m_stdSMatDataDBDB
void GetFluxDerivJacDirctnElmt(const int nConvectiveFields, const int nElmtPnt, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const Array< OneD, NekDouble > &locmu, const Array< OneD, const Array< OneD, NekDouble >> &locnormal, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
void MatrixMultiplyMatrixFreeCoeff(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out, const bool &flag=false)
virtual void v_MinusDiffusionFluxJacPoint(const int nConvectiveFields, const int nElmtPnt, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const TensorOfArray3D< NekDouble > &locDerv, const Array< OneD, NekDouble > &locmu, const Array< OneD, NekDouble > &locDmuDT, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
void CalcPreconMatBRJCoeff(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, SNekBlkMatSharedPtr >> &gmtxarray, SNekBlkMatSharedPtr &gmtVar, Array< OneD, SNekBlkMatSharedPtr > &TraceJac, Array< OneD, SNekBlkMatSharedPtr > &TraceJacDeriv, Array< OneD, Array< OneD, NekSingle >> &TraceJacDerivSign, TensorOfArray4D< NekSingle > &TraceJacArray, TensorOfArray4D< NekSingle > &TraceJacDerivArray, TensorOfArray5D< NekSingle > &TraceIPSymJacArray)
void TransTraceJacMatToArray(const Array< OneD, TypeNekBlkMatSharedPtr > &TraceJac, const Array< OneD, TypeNekBlkMatSharedPtr > &TraceJacDeriv, TensorOfArray4D< DataType > &TraceJacArray, TensorOfArray4D< DataType > &TraceJacDerivArray)
void AddMatNSBlkDiagBnd(const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield, TensorOfArray2D< TypeNekBlkMatSharedPtr > &gmtxarray, Array< OneD, TypeNekBlkMatSharedPtr > &TraceJac, Array< OneD, TypeNekBlkMatSharedPtr > &TraceJacDeriv, Array< OneD, Array< OneD, DataType >> &TraceJacDerivSign, TensorOfArray5D< DataType > &TraceIPSymJacArray)
void CalcPhysDeriv(const Array< OneD, const Array< OneD, NekDouble >> &inarray, TensorOfArray3D< NekDouble > &qfield)
void MinusDiffusionFluxJacPoint(const int nConvectiveFields, const int nElmtPnt, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const TensorOfArray3D< NekDouble > &locDerv, const Array< OneD, NekDouble > &locmu, const Array< OneD, NekDouble > &locDmuDT, const Array< OneD, NekDouble > &normals, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble >> &inpnts, Array< OneD, Array< OneD, NekDouble >> &outpnt, const NekDouble time, const NekDouble lambda)
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CFSImplicit class.
void Fill1DArrayOfBlkDiagonalMat(Array< OneD, TypeNekBlkMatSharedPtr > &gmtxarray, const DataType valu)
void CalcMuDmuDT(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &mu, Array< OneD, NekDouble > &DmuDT)
virtual void v_GetFluxDerivJacDirctnElmt(const int nConvectiveFields, const int nElmtPnt, const int nDervDir, const Array< OneD, const Array< OneD, NekDouble >> &locVars, const Array< OneD, NekDouble > &locmu, const Array< OneD, const Array< OneD, NekDouble >> &locnormal, DNekMatSharedPtr &wspMat, Array< OneD, Array< OneD, NekDouble >> &PntJacArray)
virtual void v_InitObject(bool DeclareFields=true)
Initialization object for CompressibleFlowSystem class.
SolverUtils::DiffusionSharedPtr m_diffusion
std::vector< CFSBndCondSharedPtr > m_bndConds
void GetElmtTimeStep(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, NekDouble > &tstep)
Calculate the maximum timestep on each element subject to CFL restrictions.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Compute the projection and call the method for imposing the boundary conditions in case of discontinu...
VariableConverterSharedPtr m_varConv
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
NekDouble m_NonlinIterTolRelativeL2
Definition: NekSys.h:203
NekDouble m_LinSysRelativeTolInNonlin
Definition: NekSys.h:204
void DefineNekSysResEval(FuncPointerT func, ObjectPointerT obj)
Definition: NekSys.h:100
void DefineNekSysLhsEval(FuncPointerT func, ObjectPointerT obj)
Definition: NekSys.h:107
void DefineNekSysPrecon(FuncPointerT func, ObjectPointerT obj)
Definition: NekSys.h:114
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:73
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void DefineCalcPreconMatBRJCoeff(FuncPointerT func, ObjectPointerT obj)
Definition: PreconCfsOp.h:102
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
int m_spacedim
Spatial dimension (>= expansion dim).
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
LibUtilities::CommSharedPtr m_comm
Communicator.
NekDouble m_timestep
Time step size.
NekDouble m_time
Current time of simulation.
SOLVER_UTILS_EXPORT int GetTraceTotPoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
NekDouble m_fintime
Finish time of the simulation.
NekDouble m_checktime
Time between checkpoints.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
enum HomogeneousType m_HomogeneousType
SOLVER_UTILS_EXPORT int GetNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
SOLVER_UTILS_EXPORT int GetPhys_Offset(int n)
SOLVER_UTILS_EXPORT int GetTotPoints()
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_TimeIntegLambda
coefff of spacial derivatives(rhs or m_F in GLM) in calculating the residual of the whole equation(us...
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
Array< OneD, NekDouble > m_magnitdEstimat
estimate the magnitude of each conserved varibles
static void Ger(const int &m, const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy, double *a, const int &lda)
BLAS level 2: Matrix vector multiply A = alpha*x*y**T + A where A[m x n].
Definition: Blas.hpp:211
NekNonlinSysFactory & GetNekNonlinSysFactory()
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< AssemblyMapDG > AssemblyMapDGSharedPtr
Definition: AssemblyMapDG.h:47
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< LocTraceToTraceMap > LocTraceToTraceMapSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< SNekBlkMat > SNekBlkMatSharedPtr
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
static Array< OneD, Array< OneD, Array< OneD, NekDouble > > > NullNekDoubleTensorOfArray3D
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
PreconCfsOpFactory & GetPreconCfsOpFactory()
Declaration of the boundary condition factory singleton.
Definition: PreconCfsOp.cpp:42
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:622
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
T Dot(int n, const T *w, const T *x)
vvtvp (vector times vector times vector): z = w*x*y
Definition: Vmath.cpp:1100
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.cpp:359
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.cpp:248
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:384
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
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.cpp:419
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:295
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291