CompressibleFlowSystem.cpp 33.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
///////////////////////////////////////////////////////////////////////////////
//
// File CompressibleFlowSystem.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
Douglas Serson's avatar
Douglas Serson committed
32
// Description: Compressible flow system base class with auxiliary functions
33
34
35
//
///////////////////////////////////////////////////////////////////////////////

36
#include <CompressibleFlowSolver/EquationSystems/CompressibleFlowSystem.h>
37
#include <LocalRegions/TriExp.h>
38
#include <MultiRegions/ExpList.h>
Douglas Serson's avatar
Douglas Serson committed
39

40

Chris Cantwell's avatar
Chris Cantwell committed
41
42
using namespace std;

43
44
namespace Nektar
{
45
46
47
    CompressibleFlowSystem::CompressibleFlowSystem(
        const LibUtilities::SessionReaderSharedPtr& pSession)
        : UnsteadySystem(pSession)
48
    {
49
    }
50

51
    /**
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
52
     * @brief Initialization object for CompressibleFlowSystem class.
53
     */
54
55
56
    void CompressibleFlowSystem::v_InitObject()
    {
        UnsteadySystem::v_InitObject();
57

58
59
60
        m_varConv = MemoryManager<VariableConverter>::AllocateSharedPtr(
                    m_session, m_spacedim);

61
62
        ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
                 "No UPWINDTYPE defined in session.");
63

64
65
66
        // Do not forwards transform initial condition
        m_homoInitialFwd = false;

67
        // Set up locations of velocity vector.
68
69
        m_vecLocs = Array<OneD, Array<OneD, NekDouble> >(1);
        m_vecLocs[0] = Array<OneD, NekDouble>(m_spacedim);
70
        for (int i = 0; i < m_spacedim; ++i)
71
        {
72
            m_vecLocs[0][i] = 1 + i;
73
        }
74

Douglas Serson's avatar
Douglas Serson committed
75
76
77
78
        // Loading parameters from session file
        InitialiseParameters();

        // Setting up advection and diffusion operators
79
80
81
82
83
84
85
86
87
88
89
        InitAdvection();

        // Create artificial diffusion
        if (m_shockCaptureType != "Off")
        {
            m_artificialDiffusion = GetArtificialDiffusionFactory()
                                    .CreateInstance(m_shockCaptureType,
                                                    m_session,
                                                    m_fields,
                                                    m_spacedim);
        }
Douglas Serson's avatar
Douglas Serson committed
90
91
92
93
94

        // Forcing terms for the sponge region
        m_forcing = SolverUtils::Forcing::Load(m_session, m_fields,
                                               m_fields.num_elements());

95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
        // User-defined boundary conditions
        int cnt = 0;
        for (int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
        {
            std::string type = 
                    m_fields[0]->GetBndConditions()[n]->GetUserDefined();
            if(!type.empty())
            {
                m_bndConds.push_back(GetCFSBndCondFactory().CreateInstance(
                        type,
                        m_session,
                        m_fields,
                        m_traceNormals,
                        m_spacedim,
                        n,
                        cnt));
            }
            cnt += m_fields[0]->GetBndCondExpansions()[n]->GetExpSize();
        }

Douglas Serson's avatar
Douglas Serson committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
        if (m_explicitAdvection)
        {
            m_ode.DefineOdeRhs    (&CompressibleFlowSystem::DoOdeRhs, this);
            m_ode.DefineProjection(&CompressibleFlowSystem::DoOdeProjection, this);
        }
        else
        {
            ASSERTL0(false, "Implicit CFS not set up.");
        }
    }

    /**
     * @brief Destructor for CompressibleFlowSystem class.
     */
    CompressibleFlowSystem::~CompressibleFlowSystem()
    {

    }

    /**
     * @brief Load CFS parameters from the session file
     */
    void CompressibleFlowSystem::InitialiseParameters()
    {
Douglas Serson's avatar
Douglas Serson committed
139
140
        NekDouble velInf, gasConstant;

141
142
        // Get gamma parameter from session file.
        m_session->LoadParameter("Gamma", m_gamma, 1.4);
143

Douglas Serson's avatar
Douglas Serson committed
144
145
146
147
148
        // Get gas constant from session file and compute Cp
        m_session->LoadParameter ("GasConstant",   gasConstant,   287.058);
        m_Cp      = m_gamma / (m_gamma - 1.0) * gasConstant;

        // Get pInf parameter from session file.
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
149
        m_session->LoadParameter("pInf", m_pInf, 101325);
150

Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
151
152
        // Get rhoInf parameter from session file.
        m_session->LoadParameter("rhoInf", m_rhoInf, 1.225);
153

Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
154
        // Get uInf parameter from session file.
155
        m_session->LoadParameter("uInf", velInf, 0.1);
156

157
        m_UInf = velInf*velInf;
158

Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
159
        // Get vInf parameter from session file.
160
        if (m_spacedim == 2 || m_spacedim == 3)
161
        {
162
163
            m_session->LoadParameter("vInf", velInf, 0.0);
            m_UInf += velInf*velInf;
164
        }
165

Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
166
        // Get wInf parameter from session file.
167
        if (m_spacedim == 3)
168
        {
169
170
            m_session->LoadParameter("wInf", velInf, 0.0);
            m_UInf += velInf*velInf;
171
        }
172
        m_UInf = sqrt(m_UInf);
173

Douglas Serson's avatar
Douglas Serson committed
174
        // Viscosity
175
176
        m_session->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
        m_session->LoadParameter ("mu",            m_mu,            1.78e-05);
Douglas Serson's avatar
Douglas Serson committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195

        // Thermal conductivity or Prandtl
        if( m_session->DefinesParameter("thermalConductivity"))
        {
            ASSERTL0( !m_session->DefinesParameter("Pr"),
                 "Cannot define both Pr and thermalConductivity.");

            m_session->LoadParameter ("thermalConductivity",
                                        m_thermalConductivity);
            m_Prandtl = m_Cp * m_mu / m_thermalConductivity;
        }
        else
        {
            m_session->LoadParameter ("Pr",
                                        m_Prandtl, 0.72);
            m_thermalConductivity = m_Cp * m_mu / m_Prandtl;
        }

        // Steady state tolerance
196
        m_session->LoadParameter("SteadyStateTol", m_steadyStateTol, 0.0);
197
198
199
200

        // Shock capture
        m_session->LoadSolverInfo("ShockCaptureType",
                                  m_shockCaptureType,    "Off");
Douglas Serson's avatar
Douglas Serson committed
201
    }
Dave Moxey's avatar
Dave Moxey committed
202

Douglas Serson's avatar
Douglas Serson committed
203
204
205
    /**
     * @brief Create advection and diffusion objects for CFS
     */
206
    void CompressibleFlowSystem::InitAdvection()
Douglas Serson's avatar
Douglas Serson committed
207
208
209
210
    {
        // Check if projection type is correct
        ASSERTL0(m_projectionType == MultiRegions::eDiscontinuous,
                "Unsupported projection type.");
211

212
        string advName, riemName;
Douglas Serson's avatar
Douglas Serson committed
213
        m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
214

Douglas Serson's avatar
Douglas Serson committed
215
216
217
218
        m_advection = SolverUtils::GetAdvectionFactory()
                                    .CreateInstance(advName, advName);

        if (m_specHP_dealiasing)
219
        {
Douglas Serson's avatar
Douglas Serson committed
220
221
            m_advection->SetFluxVector(&CompressibleFlowSystem::
                                       GetFluxVectorDeAlias, this);
222
223
224
        }
        else
        {
Douglas Serson's avatar
Douglas Serson committed
225
226
            m_advection->SetFluxVector  (&CompressibleFlowSystem::
                                          GetFluxVector, this);
227
        }
228

Douglas Serson's avatar
Douglas Serson committed
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
        // Setting up Riemann solver for advection operator
        m_session->LoadSolverInfo("UpwindType", riemName, "Average");

        SolverUtils::RiemannSolverSharedPtr riemannSolver;
        riemannSolver = SolverUtils::GetRiemannSolverFactory()
                                    .CreateInstance(riemName);

        // Setting up parameters for advection operator Riemann solver
        riemannSolver->SetParam (
            "gamma",   &CompressibleFlowSystem::GetGamma,   this);
        riemannSolver->SetAuxVec(
            "vecLocs", &CompressibleFlowSystem::GetVecLocs, this);
        riemannSolver->SetVector(
            "N",       &CompressibleFlowSystem::GetNormals, this);

        // Concluding initialisation of advection / diffusion operators
        m_advection->SetRiemannSolver   (riemannSolver);
        m_advection->InitObject         (m_session, m_fields);
247
    }
248

249
250
251
252
253
254
255
256
257
258
259
    /**
     * @brief Compute the right-hand side.
     */
    void CompressibleFlowSystem::DoOdeRhs(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray,
              Array<OneD,       Array<OneD, NekDouble> > &outarray,
        const NekDouble                                   time)
    {
        int i;
        int nvariables = inarray.num_elements();
        int npoints    = GetNpoints();
260
        int nTracePts  = GetTraceTotPoints();
261

262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
        // Store forwards/backwards space along trace space
        Array<OneD, Array<OneD, NekDouble> > Fwd    (nvariables);
        Array<OneD, Array<OneD, NekDouble> > Bwd    (nvariables);

        if (m_HomogeneousType == eHomogeneous1D)
        {
            Fwd = NullNekDoubleArrayofArray;
            Bwd = NullNekDoubleArrayofArray;
        }
        else
        {
            for(i = 0; i < nvariables; ++i)
            {
                Fwd[i]     = Array<OneD, NekDouble>(nTracePts, 0.0);
                Bwd[i]     = Array<OneD, NekDouble>(nTracePts, 0.0);
                m_fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
            }
        }

        // Calculate advection
        DoAdvection(inarray, outarray, time, Fwd, Bwd);
283
284
285
286
287
288
289

        // Negate results
        for (i = 0; i < nvariables; ++i)
        {
            Vmath::Neg(npoints, outarray[i], 1);
        }

290
291
        // Add diffusion terms
        DoDiffusion(inarray, outarray, Fwd, Bwd);
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345

        // Add forcing terms
        std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
        for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
        {
            (*x)->Apply(m_fields, inarray, outarray, time);
        }
    }

    /**
     * @brief Compute the projection and call the method for imposing the 
     * boundary conditions in case of discontinuous projection.
     */
    void CompressibleFlowSystem::DoOdeProjection(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray,
              Array<OneD,       Array<OneD, NekDouble> > &outarray,
        const NekDouble                                   time)
    {
        int i;
        int nvariables = inarray.num_elements();

        switch(m_projectionType)
        {
            case MultiRegions::eDiscontinuous:
            {
                // Just copy over array
                int npoints = GetNpoints();

                for(i = 0; i < nvariables; ++i)
                {
                    Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
                }
                SetBoundaryConditions(outarray, time);
                break;
            }
            case MultiRegions::eGalerkin:
            case MultiRegions::eMixed_CG_Discontinuous:
            {
                ASSERTL0(false, "No Continuous Galerkin for full compressible "
                                "Navier-Stokes equations");
                break;
            }
            default:
                ASSERTL0(false, "Unknown projection scheme");
                break;
        }
    }

    /**
     * @brief Compute the advection terms for the right-hand side
     */
    void CompressibleFlowSystem::DoAdvection(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray,
              Array<OneD,       Array<OneD, NekDouble> > &outarray,
346
347
348
        const NekDouble                                   time,
        const Array<OneD, Array<OneD, NekDouble> >       &pFwd,
        const Array<OneD, Array<OneD, NekDouble> >       &pBwd)
349
350
351
352
353
    {
        int nvariables = inarray.num_elements();
        Array<OneD, Array<OneD, NekDouble> > advVel(m_spacedim);

        m_advection->Advect(nvariables, m_fields, advVel, inarray,
354
                            outarray, time, pFwd, pBwd);
355
356
357
358
359
360
361
    }

    /**
     * @brief Add the diffusions terms to the right-hand side
     */
    void CompressibleFlowSystem::DoDiffusion(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray,
362
363
364
              Array<OneD,       Array<OneD, NekDouble> > &outarray,
            const Array<OneD, Array<OneD, NekDouble> >   &pFwd,
            const Array<OneD, Array<OneD, NekDouble> >   &pBwd)
365
    {
366
        v_DoDiffusion(inarray, outarray, pFwd, pBwd);
367
368
369
370
371

        if (m_shockCaptureType != "Off")
        {
            m_artificialDiffusion->DoArtificialDiffusion(inarray, outarray);
        }
372
373
    }

374
375
376
377
    void CompressibleFlowSystem::SetBoundaryConditions(
            Array<OneD, Array<OneD, NekDouble> >             &physarray,
            NekDouble                                         time)
    {
378
379
        int nTracePts  = GetTraceTotPoints();
        int nvariables = physarray.num_elements();
380

381
382
383
384
385
386
        Array<OneD, Array<OneD, NekDouble> > Fwd(nvariables);
        for (int i = 0; i < nvariables; ++i)
        {
            Fwd[i] = Array<OneD, NekDouble>(nTracePts);
            m_fields[i]->ExtractTracePhys(physarray[i], Fwd[i]);
        }
387

388
389
        if (m_bndConds.size())
        {
390
391
392
393
394
395
            // Loop over user-defined boundary conditions
            std::vector<CFSBndCondSharedPtr>::iterator x;
            for (x = m_bndConds.begin(); x != m_bndConds.end(); ++x)
            {
                (*x)->Apply(Fwd, physarray, time);
            }
396
397
        }
    }
398

399
400
    /**
     * @brief Return the flux vector for the compressible Euler equations.
401
     *
402
403
404
405
     * @param physfield   Fields.
     * @param flux        Resulting flux.
     */
    void CompressibleFlowSystem::GetFluxVector(
Dave Moxey's avatar
Dave Moxey committed
406
407
        const Array<OneD, Array<OneD, NekDouble> >               &physfield,
              Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
408
    {
409
        int i, j;
410
        int nq = physfield[0].num_elements();
411
        int nVariables = m_fields.num_elements();
412

413
414
        Array<OneD, NekDouble> pressure(nq);
        Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
415

416
417
418
419
420
421
        // Flux vector for the rho equation
        for (i = 0; i < m_spacedim; ++i)
        {
            velocity[i] = Array<OneD, NekDouble>(nq);
            Vmath::Vcopy(nq, physfield[i+1], 1, flux[0][i], 1);
        }
422

423
424
        m_varConv->GetVelocityVector(physfield, velocity);
        m_varConv->GetPressure(physfield, velocity, pressure);
425

426
427
428
429
        // Flux vector for the velocity fields
        for (i = 0; i < m_spacedim; ++i)
        {
            for (j = 0; j < m_spacedim; ++j)
430
            {
431
432
                Vmath::Vmul(nq, velocity[j], 1, physfield[i+1], 1,
                            flux[i+1][j], 1);
433
            }
434

435
436
437
            // Add pressure to appropriate field
            Vmath::Vadd(nq, flux[i+1][i], 1, pressure, 1, flux[i+1][i], 1);
        }
438

439
        // Flux vector for energy.
440
        Vmath::Vadd(nq, physfield[m_spacedim+1], 1, pressure, 1,
441
                    pressure, 1);
442

443
444
        for (j = 0; j < m_spacedim; ++j)
        {
445
            Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
446
447
                        flux[m_spacedim+1][j], 1);
        }
448

449
450
        // For the smooth viscosity model
        if (nVariables == m_spacedim+3)
451
        {
Dave Moxey's avatar
Dave Moxey committed
452
            // Add a zero row for the advective fluxes
453
454
            for (j = 0; j < m_spacedim; ++j)
            {
Dave Moxey's avatar
Dave Moxey committed
455
                Vmath::Zero(nq, flux[m_spacedim+2][j], 1);
456
457
458
            }
        }
    }
459

460
461
462
    /**
     * @brief Return the flux vector for the compressible Euler equations
     * by using the de-aliasing technique.
463
     *
464
465
466
467
468
469
470
471
     * @param physfield   Fields.
     * @param flux        Resulting flux.
     */
    void CompressibleFlowSystem::GetFluxVectorDeAlias(
        const Array<OneD, Array<OneD, NekDouble> >               &physfield,
              Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux)
    {
        int i, j;
472
        int nq = physfield[0].num_elements();
473
        int nVariables = m_fields.num_elements();
474

475
        // Factor to rescale 1d points in dealiasing
476
477
        NekDouble OneDptscale = 2;
        nq = m_fields[0]->Get1DScaledTotPoints(OneDptscale);
478

479
480
        Array<OneD, NekDouble> pressure(nq);
        Array<OneD, Array<OneD, NekDouble> > velocity(m_spacedim);
481

482
        Array<OneD, Array<OneD, NekDouble> > physfield_interp(nVariables);
483
        Array<OneD, Array<OneD, Array<OneD, NekDouble> > > flux_interp(
484
                                                            nVariables);
485

486
        for (i = 0; i < nVariables; ++ i)
487
488
489
        {
            physfield_interp[i] = Array<OneD, NekDouble>(nq);
            flux_interp[i] = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
490
491
            m_fields[0]->PhysInterp1DScaled(
                OneDptscale, physfield[i], physfield_interp[i]);
492

493
            for (j = 0; j < m_spacedim; ++j)
494
            {
495
                flux_interp[i][j] = Array<OneD, NekDouble>(nq);
496
            }
497
        }
498

499
500
501
502
        // Flux vector for the rho equation
        for (i = 0; i < m_spacedim; ++i)
        {
            velocity[i] = Array<OneD, NekDouble>(nq);
503

504
            // Galerkin project solution back to original space
505
506
            m_fields[0]->PhysGalerkinProjection1DScaled(
                OneDptscale, physfield_interp[i+1], flux[0][i]);
507
        }
508

509
510
        m_varConv->GetVelocityVector(physfield_interp, velocity);
        m_varConv->GetPressure      (physfield_interp, velocity, pressure);
511

512
513
514
        // Evaluation of flux vector for the velocity fields
        for (i = 0; i < m_spacedim; ++i)
        {
515
            for (j = 0; j < m_spacedim; ++j)
516
            {
517
518
519
                Vmath::Vmul(nq, velocity[j], 1, physfield_interp[i+1], 1,
                            flux_interp[i+1][j], 1);
            }
520

521
            // Add pressure to appropriate field
522
            Vmath::Vadd(nq, flux_interp[i+1][i], 1, pressure,1,
523
524
                        flux_interp[i+1][i], 1);
        }
525

526
527
528
529
530
        // Galerkin project solution back to origianl space
        for (i = 0; i < m_spacedim; ++i)
        {
            for (j = 0; j < m_spacedim; ++j)
            {
531
532
                m_fields[0]->PhysGalerkinProjection1DScaled(
                    OneDptscale, flux_interp[i+1][j], flux[i+1][j]);
533
534
            }
        }
535

536
        // Evaluation of flux vector for energy
537
        Vmath::Vadd(nq, physfield_interp[m_spacedim+1], 1, pressure, 1,
538
                    pressure, 1);
539

540
541
        for (j = 0; j < m_spacedim; ++j)
        {
542
            Vmath::Vmul(nq, velocity[j], 1, pressure, 1,
543
                        flux_interp[m_spacedim+1][j], 1);
544

545
            // Galerkin project solution back to origianl space
546
547
548
549
            m_fields[0]->PhysGalerkinProjection1DScaled(
                OneDptscale,
                flux_interp[m_spacedim+1][j],
                flux[m_spacedim+1][j]);
550
551
552
        }
    }

553
    /**
Dave Moxey's avatar
Dave Moxey committed
554
555
     * @brief Perform post-integration checks, presently just to check steady
     * state behaviour.
556
     */
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
557
    bool CompressibleFlowSystem::v_PostIntegrate(int step)
558
    {
559
        if (m_steadyStateTol > 0.0)
560
        {
561
562
            bool doOutput = step % m_infosteps == 0;
            if (CalcSteadyState(doOutput))
563
564
565
566
567
568
569
570
571
            {
                if (m_comm->GetRank() == 0)
                {
                    cout << "Reached Steady State to tolerance "
                         << m_steadyStateTol << endl;
                }
                return true;
            }
        }
572
573
        return false;
    }
Dave Moxey's avatar
Dave Moxey committed
574
575
576
577
578

    /**
     * @brief Calculate whether the system has reached a steady state by
     * observing residuals to a user-defined tolerance.
     */
579
    bool CompressibleFlowSystem::CalcSteadyState(bool output)
580
    {
581
582
583
584
        const int nPoints = GetTotPoints();
        const int nFields = m_fields.num_elements();

        // Holds L2 errors.
Dave Moxey's avatar
Dave Moxey committed
585
586
        Array<OneD, NekDouble> L2      (nFields);
        Array<OneD, NekDouble> residual(nFields);
587
588
589
590
591
592
593
594

        for (int i = 0; i < nFields; ++i)
        {
            Array<OneD, NekDouble> diff(nPoints);

            Vmath::Vsub(nPoints, m_fields[i]->GetPhys(), 1, m_un[i], 1, diff, 1);
            Vmath::Vmul(nPoints, diff, 1, diff, 1, diff, 1);
            residual[i] = Vmath::Vsum(nPoints, diff, 1);
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
595
596
        }

597
        m_comm->AllReduce(residual, LibUtilities::ReduceSum);
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
598

599
600
601
602
        // L2 error
        L2[0] = sqrt(residual[0]) / m_rhoInf;

        for (int i = 1; i < nFields-1; ++i)
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
603
        {
604
            L2[i] = sqrt(residual[i]) / m_UInf / m_rhoInf;
605
        }
606
607
608
609
610

        NekDouble Einf = m_pInf / (m_gamma-1.0) + 0.5 * m_rhoInf * m_UInf;
        L2[nFields-1] = sqrt(residual[nFields-1]) / Einf;

        if (m_comm->GetRank() == 0 && output)
611
        {
612
613
614
615
616
            // Output time
            m_errFile << setprecision(8) << setw(17) << scientific << m_time;

            // Output residuals
            for (int i = 0; i < nFields; ++i)
617
            {
618
619
                m_errFile << setprecision(11) << setw(22) << scientific
                          << L2[i];
620
            }
621
622

            m_errFile << endl;
623
        }
624
625
626
627
628
629

        // Calculate maximum L2 error
        NekDouble maxL2 = Vmath::Vmax(nFields, L2, 1);

        if (m_session->DefinesCmdLineArgument("verbose") &&
            m_comm->GetRank() == 0 && output)
630
        {
631
            cout << "-- Maximum L^2 residual: " << maxL2 << endl;
632
        }
633

634
635
        if (maxL2 <= m_steadyStateTol)
        {
636
            return true;
637
        }
638

Dave Moxey's avatar
Dave Moxey committed
639
640
641
642
643
        for (int i = 0; i < m_fields.num_elements(); ++i)
        {
            Vmath::Vcopy(nPoints, m_fields[i]->GetPhys(), 1, m_un[i], 1);
        }

644
        return false;
645
    }
646

647
648
649
    /**
     * @brief Calculate the maximum timestep subject to CFL restrictions.
     */
650
651
    NekDouble CompressibleFlowSystem::v_GetTimeStep(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray)
652
    {
653
        int n;
654
        int nElements = m_fields[0]->GetExpSize();
655
        const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();
656

657
658
        Array<OneD, NekDouble> tstep      (nElements, 0.0);
        Array<OneD, NekDouble> stdVelocity(nElements);
659

660
        // Get standard velocity to compute the time-step limit
661
        GetStdVelocity(inarray, stdVelocity);
662

663
        // Factors to compute the time-step limit
Chris Cantwell's avatar
Chris Cantwell committed
664
665
666
        NekDouble minLength = 0.0;
        NekDouble alpha     = MaxTimeStepEstimator();
        NekDouble cLambda   = 0.2; // Spencer book-317
667

668
669
670
671
672
673
        // Loop over elements to compute the time-step limit for each element
        for(n = 0; n < nElements; ++n)
        {
            int npoints = m_fields[0]->GetExp(n)->GetTotPoints();
            Array<OneD, NekDouble> one2D(npoints, 1.0);
            NekDouble Area = m_fields[0]->GetExp(n)->Integral(one2D);
674

Chris Cantwell's avatar
Chris Cantwell committed
675
            minLength = sqrt(Area);
676
            if (m_fields[0]->GetExp(n)->as<LocalRegions::TriExp>())
677
            {
Chris Cantwell's avatar
Chris Cantwell committed
678
                minLength *= 2.0;
Dirk Ekelschot's avatar
Dirk Ekelschot committed
679
680
            }

681
682
            tstep[n] = m_cflSafetyFactor * alpha * minLength
                     / (stdVelocity[n] * cLambda
683
                        * (ExpOrder[n] - 1) * (ExpOrder[n] - 1));
684
        }
685

686
687
688
689
690
        // Get the minimum time-step limit and return the time-step
        NekDouble TimeStep = Vmath::Vmin(nElements, tstep, 1);
        m_comm->AllReduce(TimeStep, LibUtilities::ReduceMin);
        return TimeStep;
    }
691

692
693
694
695
696
697
698
    /**
     * @brief Set up logic for residual calculation.
     */
    void CompressibleFlowSystem::v_SetInitialConditions(
        NekDouble initialtime,
        bool      dumpInitialConditions,
        const int domain)
699
700
701
702
703
704
705
706
707
708
709
710
711
    {
        EquationSystem::v_SetInitialConditions(initialtime, false);

        // insert white noise in initial condition
        NekDouble Noise;
        int phystot = m_fields[0]->GetTotPoints();
        Array<OneD, NekDouble> noise(phystot);

        m_session->LoadParameter("Noise", Noise,0.0);
        int m_nConvectiveFields =  m_fields.num_elements();

        if (Noise > 0.0)
        {
Douglas Serson's avatar
Douglas Serson committed
712
            int seed = - (m_comm->GetRank() * m_nConvectiveFields);
713
714
715
            for (int i = 0; i < m_nConvectiveFields; i++)
            {
                Vmath::FillWhiteNoise(phystot, Noise, noise, 1,
Douglas Serson's avatar
Douglas Serson committed
716
717
718
                                      seed);
                --seed;

719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
                Vmath::Vadd(phystot, m_fields[i]->GetPhys(), 1,
                            noise, 1, m_fields[i]->UpdatePhys(), 1);
                m_fields[i]->FwdTrans_IterPerExp(m_fields[i]->GetPhys(),
                                                 m_fields[i]->UpdateCoeffs());
            }
        }

        InitializeSteadyState();

        if (dumpInitialConditions)
        {
            // Dump initial conditions to file
            Checkpoint_Output(0);
        }
    }

    void CompressibleFlowSystem::InitializeSteadyState()
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
    {
        if (m_session->DefinesParameter("SteadyStateTol"))
        {
            const int nPoints = m_fields[0]->GetTotPoints();
            m_un = Array<OneD, Array<OneD, NekDouble> > (
                m_fields.num_elements());

            for (int i = 0; i < m_fields.num_elements(); ++i)
            {
                m_un[i] = Array<OneD, NekDouble>(nPoints);
                Vmath::Vcopy(nPoints, m_fields[i]->GetPhys(), 1, m_un[i], 1);
            }

            if (m_comm->GetRank() == 0)
            {
                std::string fName = m_session->GetSessionName() +
                    std::string(".res");
                m_errFile.open(fName.c_str());
                m_errFile << "# "
                          << setw(15) << left << "Time"
                          << setw(22) << left << "rho";

                std::string velFields[3] = {"u", "v", "w"};

                for (int i = 0; i < m_fields.num_elements()-2; ++i)
                {
                    m_errFile << setw(22) << "rho"+velFields[i];
                }

                m_errFile << setw(22) << left << "E" << endl;
            }
        }
    }


771
    /**
772
     * @brief Compute the advection velocity in the standard space
773
     * for each element of the expansion.
774
     *
775
776
777
     * @param inarray    Momentum field.
     * @param stdV       Standard velocity field.
     */
778
779
780
    void CompressibleFlowSystem::GetStdVelocity(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray,
              Array<OneD,                   NekDouble>   &stdV)
781
    {
782
783
        int nTotQuadPoints = GetTotPoints();
        int n_element      = m_fields[0]->GetExpSize();
784
        int nBCEdgePts           = 0;
785
786

        // Getting the velocity vector on the 2D normal space
787
788
        Array<OneD, Array<OneD, NekDouble> > velocity   (m_spacedim);
        Array<OneD, Array<OneD, NekDouble> > stdVelocity(m_spacedim);
789
790
        Array<OneD, NekDouble>               pressure   (nTotQuadPoints);
        Array<OneD, NekDouble>               soundspeed (nTotQuadPoints);
Chris Cantwell's avatar
Chris Cantwell committed
791
        LibUtilities::PointsKeyVector        ptsKeys;
792

793
        // Zero output array
794
        Vmath::Zero(stdV.num_elements(), stdV, 1);
795

796
        for (int i = 0; i < m_spacedim; ++i)
797
798
799
800
        {
            velocity   [i] = Array<OneD, NekDouble>(nTotQuadPoints);
            stdVelocity[i] = Array<OneD, NekDouble>(nTotQuadPoints, 0.0);
        }
Dave Moxey's avatar
Dave Moxey committed
801

802
803
804
        m_varConv->GetVelocityVector(inarray, velocity);
        m_varConv->GetPressure      (inarray, velocity, pressure);
        m_varConv->GetSoundSpeed    (inarray, pressure, soundspeed);
805

806
        for(int el = 0; el < n_element; ++el)
807
        {
Chris Cantwell's avatar
Chris Cantwell committed
808
809
            ptsKeys = m_fields[0]->GetExp(el)->GetPointsKeys();

810
            // Possible bug: not multiply by jacobian??
811
812
            const SpatialDomains::GeomFactorsSharedPtr metricInfo =
                m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo();
813
            const Array<TwoD, const NekDouble> &gmat =
814
815
                m_fields[0]->GetExp(el)->GetGeom()->GetMetricInfo()
                                                  ->GetDerivFactors(ptsKeys);
816

817
            int nq = m_fields[0]->GetExp(el)->GetTotPoints();
818

819
            if(metricInfo->GetGtype() == SpatialDomains::eDeformed)
820
821
            {
                // d xi/ dx = gmat = 1/J * d x/d xi
822
                for (int i = 0; i < m_spacedim; ++i)
823
                {
824
825
826
                    Vmath::Vmul(nq, gmat[i], 1, velocity[0], 1,
                                stdVelocity[i], 1);
                    for (int j = 1; j < m_spacedim; ++j)
827
                    {
828
                        Vmath::Vvtvp(nq, gmat[m_spacedim*j+i], 1, velocity[j],
829
                                     1, stdVelocity[i], 1, stdVelocity[i], 1);
830
831
832
833
834
                    }
                }
            }
            else
            {
835
                for (int i = 0; i < m_spacedim; ++i)
836
                {
837
838
839
                    Vmath::Smul(nq, gmat[i][0], velocity[0], 1,
                                stdVelocity[i], 1);
                    for (int j = 1; j < m_spacedim; ++j)
840
                    {
841
                        Vmath::Svtvp(nq, gmat[m_spacedim*j+i][0], velocity[j],
842
                                     1, stdVelocity[i], 1, stdVelocity[i], 1);
843
844
845
                    }
                }
            }
846

847
            for (int i = 0; i < nq; ++i)
848
            {
849
                NekDouble pntVelocity = 0.0;
850
                for (int j = 0; j < m_spacedim; ++j)
851
852
853
                {
                    pntVelocity += stdVelocity[j][i]*stdVelocity[j][i];
                }
854
                pntVelocity = sqrt(pntVelocity) + soundspeed[nBCEdgePts];
855
856
                if (pntVelocity > stdV[el])
                {
857
                    stdV[el] = pntVelocity;
858
                }
859
                nBCEdgePts++;
860
861
862
            }
        }
    }
863

864
    /**
865
866
     * @brief Set the denominator to compute the time step when a cfl
     * control is employed. This function is no longer used but is still
867
868
869
870
871
872
     * here for being utilised in the future.
     *
     * @param n   Order of expansion element by element.
     */
    NekDouble CompressibleFlowSystem::GetStabilityLimit(int n)
    {
873
874
875
        ASSERTL0(n <= 20, "Illegal modes dimension for CFL calculation "
                          "(P has to be less then 20)");

876
877
        NekDouble CFLDG[21] = {  2.0000,   6.0000,  11.8424,  19.1569,
                                27.8419,  37.8247,  49.0518,  61.4815,
878
879
880
881
                                75.0797,  89.8181, 105.6700, 122.6200,
                               140.6400, 159.7300, 179.8500, 201.0100,
                               223.1800, 246.3600, 270.5300, 295.6900,
                               321.8300}; //CFLDG 1D [0-20]
Dave Moxey's avatar
Dave Moxey committed
882
        NekDouble CFL = 0.0;
883

884
885
886
887
        if (m_projectionType == MultiRegions::eDiscontinuous)
        {
            CFL = CFLDG[n];
        }
888
        else
889
        {
890
891
            ASSERTL0(false, "Continuous Galerkin stability coefficients "
                            "not introduced yet.");
892
        }
893

894
895
        return CFL;
    }
896

897
    /**
898
899
     * @brief Compute the vector of denominators to compute the time step
     * when a cfl control is employed. This function is no longer used but
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
     * is still here for being utilised in the future.
     *
     * @param ExpOrder   Order of expansion element by element.
     */
    Array<OneD, NekDouble> CompressibleFlowSystem::GetStabilityLimitVector(
        const Array<OneD,int> &ExpOrder)
    {
        int i;
        Array<OneD,NekDouble> returnval(m_fields[0]->GetExpSize(), 0.0);
        for (i =0; i<m_fields[0]->GetExpSize(); i++)
        {
            returnval[i] = GetStabilityLimit(ExpOrder[i]);
        }
        return returnval;
    }
Dirk Ekelschot's avatar
Dirk Ekelschot committed
915

916
917
918
919
    void CompressibleFlowSystem::v_ExtraFldOutput(
        std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
        std::vector<std::string>             &variables)
    {
920
921
922
923
        bool extraFields;
        m_session->MatchSolverInfo("OutputExtraFields","True",
                                   extraFields, true);
        if (extraFields)
924
        {
925
926
927
928
929
930
931
932
933
934
            const int nPhys   = m_fields[0]->GetNpoints();
            const int nCoeffs = m_fields[0]->GetNcoeffs();
            Array<OneD, Array<OneD, NekDouble> > tmp(m_fields.num_elements());

            for (int i = 0; i < m_fields.num_elements(); ++i)
            {
                tmp[i] = m_fields[i]->GetPhys();
            }

            Array<OneD, NekDouble> pressure(nPhys), soundspeed(nPhys), mach(nPhys);
935
            Array<OneD, NekDouble> sensor(nPhys), SensorKappa(nPhys);
936

937
938
939
            m_varConv->GetPressure  (tmp, pressure);
            m_varConv->GetSoundSpeed(tmp, pressure, soundspeed);
            m_varConv->GetMach      (tmp, soundspeed, mach);
940
            m_varConv->GetSensor    (m_fields[0], tmp, sensor, SensorKappa);
941
942

            Array<OneD, NekDouble> pFwd(nCoeffs), sFwd(nCoeffs), mFwd(nCoeffs);
943
            Array<OneD, NekDouble> sensFwd(nCoeffs);
944

945
946
947
948
            m_fields[0]->FwdTrans_IterPerExp(pressure,   pFwd);
            m_fields[0]->FwdTrans_IterPerExp(soundspeed, sFwd);
            m_fields[0]->FwdTrans_IterPerExp(mach,       mFwd);
            m_fields[0]->FwdTrans_IterPerExp(sensor,     sensFwd);
949
950
951
952
953
954
955
956
957

            variables.push_back  ("p");
            variables.push_back  ("a");
            variables.push_back  ("Mach");
            variables.push_back  ("Sensor");
            fieldcoeffs.push_back(pFwd);
            fieldcoeffs.push_back(sFwd);
            fieldcoeffs.push_back(mFwd);
            fieldcoeffs.push_back(sensFwd);
958
959
960
961
962

            if (m_artificialDiffusion)
            {
                // reuse pressure
                m_artificialDiffusion->GetArtificialViscosity(tmp, pressure);
963
                m_fields[0]->FwdTrans_IterPerExp(pressure,   pFwd);
964
965
966
967

                variables.push_back  ("ArtificialVisc");
                fieldcoeffs.push_back(pFwd);
            }
968
        }
969
    }
970
}