UnsteadySystem.cpp 38 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
32
33
34
35
36
37
38
///////////////////////////////////////////////////////////////////////////////
//
// File UnsteadySystem.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.
//
// Description: Generic timestepping for Unsteady solvers
//
///////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <iomanip>

39
#include <LibUtilities/TimeIntegration/TimeIntegrationWrapper.h>
40
#include <LibUtilities/BasicUtils/Timer.h>
41
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
42
43
#include <SolverUtils/UnsteadySystem.h>

44
45
using namespace std;

46
47
48
49
50
51
52
53
namespace Nektar
{	
    namespace SolverUtils
    {
        /**
         * @class UnsteadySystem
         *
         * Provides the underlying timestepping framework for unsteady solvers
54
55
56
         * including the general timestepping routines. This class is not 
         * intended to be directly instantiated, but rather is a base class 
         * on which to define unsteady solvers.
57
58
59
60
61
62
63
64
65
66
67
68
         *
         * For details on implementing unsteady solvers see
         * \ref sectionADRSolverModuleImplementation here
         */

        /**
         * Processes SolverInfo parameters from the session file and sets up
         * timestepping-specific code.
         * @param   pSession        Session object to read parameters from.
         */
        UnsteadySystem::UnsteadySystem(
            const LibUtilities::SessionReaderSharedPtr& pSession)
69
70
71
            : EquationSystem(pSession),
              m_infosteps(10)

72
73
74
        {
        }

Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
75
        /**
76
         * Initialization object for UnsteadySystem class.
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
77
         */
78
79
80
        void UnsteadySystem::v_InitObject()
        {
            EquationSystem::v_InitObject();
81
82
            
            m_initialStep = 0;
83
84
85
86
87
88
89
90
91

            // Load SolverInfo parameters
            m_session->MatchSolverInfo("DIFFUSIONADVANCEMENT","Explicit",
                                       m_explicitDiffusion,true);
            m_session->MatchSolverInfo("ADVECTIONADVANCEMENT","Explicit",
                                       m_explicitAdvection,true);
            m_session->MatchSolverInfo("REACTIONADVANCEMENT", "Explicit",
                                       m_explicitReaction, true);

92
93
            // For steady problems, we do not initialise the time integration
            if (m_session->DefinesSolverInfo("TIMEINTEGRATIONMETHOD"))
94
            {
95
                m_intScheme = LibUtilities::GetTimeIntegrationWrapperFactory().
96
97
                    CreateInstance(m_session->GetSolverInfo(
                                       "TIMEINTEGRATIONMETHOD"));
98

99
100
101
                // Load generic input parameters
                m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
                m_session->LoadParameter("CFL", m_cflSafetyFactor, 0.0);
102

103
104
105
106
                // Set up time to be dumped in field information
                m_fieldMetaDataMap["Time"] =
                        boost::lexical_cast<std::string>(m_time);
            }
107

108
109
110
            // By default attempt to forward transform initial condition.
            m_homoInitialFwd = true;

111
112
113
114
115
            // Set up filters
            LibUtilities::FilterMap::const_iterator x;
            LibUtilities::FilterMap f = m_session->GetFilters();
            for (x = f.begin(); x != f.end(); ++x)
            {
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
116
117
118
119
                m_filters.push_back(GetFilterFactory().CreateInstance(
                                                                x->first, 
                                                                m_session, 
                                                                x->second));
120
            }
121
122
        }
        
123
        /**
124
         * Destructor for the class UnsteadyAdvection.
125
126
127
128
129
         */
        UnsteadySystem::~UnsteadySystem()
        {
        }

Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
130
        /**
131
         * @brief Returns the maximum time estimator for CFL control.
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
132
133
134
         */
        NekDouble UnsteadySystem::MaxTimeStepEstimator()
        {
135
            NekDouble TimeStability = 0.0;
136
            switch(m_intScheme->GetIntegrationMethod())
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
137
138
139
            {
                case LibUtilities::eForwardEuler:
                case LibUtilities::eClassicalRungeKutta4:
140
                case LibUtilities::eRungeKutta4:
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
141
142
143
144
145
                {
                    TimeStability = 2.784;
                    break;
                }
                case LibUtilities::eAdamsBashforthOrder1:
146
147
                case LibUtilities::eMidpoint:
                case LibUtilities::eRungeKutta2:
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
148
                case LibUtilities::eRungeKutta2_ImprovedEuler:
149
150
                case LibUtilities::eRungeKutta2_SSP:
                case LibUtilities::eRungeKutta3_SSP:
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
                {
                    TimeStability = 2.0;
                    break;
                }
                case LibUtilities::eAdamsBashforthOrder2:
                {
                    TimeStability = 1.0;
                    break;
                }
                default:
                {
                    ASSERTL0(
                        false,
                        "No CFL control implementation for this time"
                        "integration scheme");
                }
            }
            return TimeStability;
        }
        
171
        /**
172
173
         * @brief Initialises the time integration scheme (as specified in the 
         * session file), and perform the time integration.
174
175
176
         */
        void UnsteadySystem::v_DoSolve()
        {
177
178
            ASSERTL0(m_intScheme != 0, "No time integration scheme.");

179
            int i = 1;
180
            int nvariables = 0;
181
            int nfields = m_fields.num_elements();
182
183
184

            if (m_intVariables.empty())
            {
185
                for (i = 0; i < nfields; ++i)
186
187
188
                {
                    m_intVariables.push_back(i);
                }
189
                nvariables = nfields;
190
191
192
193
194
195
            }
            else
            {
                nvariables = m_intVariables.size();
            }

Chris Cantwell's avatar
Chris Cantwell committed
196
            // Integrate in wave-space if using homogeneous1D
197
            if(m_HomogeneousType == eHomogeneous1D && m_homoInitialFwd)
198
199
200
            {
                for(i = 0; i < nfields; ++i)
                {
Chris Cantwell's avatar
Chris Cantwell committed
201
202
                    m_fields[i]->HomogeneousFwdTrans(m_fields[i]->GetPhys(),
                                                     m_fields[i]->UpdatePhys());
203
204
205
206
207
                    m_fields[i]->SetWaveSpace(true);
                    m_fields[i]->SetPhysState(false);
                }
            }

208
            // Set up wrapper to fields data storage.
209
210
211
            Array<OneD, Array<OneD, NekDouble> > fields(nvariables);
            Array<OneD, Array<OneD, NekDouble> > tmp   (nvariables);
            
212
            // Order storage to list time-integrated fields first.
213
214
            for(i = 0; i < nvariables; ++i)
            {
215
                fields[i] = m_fields[m_intVariables[i]]->GetPhys();
216
217
                m_fields[m_intVariables[i]]->SetPhysState(false);
            }
218
            
Chris Cantwell's avatar
Chris Cantwell committed
219
            // Initialise time integration scheme
220
221
            m_intSoln = m_intScheme->InitializeScheme(
                m_timestep, fields, m_time, m_ode);
222

Chris Cantwell's avatar
Chris Cantwell committed
223
            // Initialise filters
224
225
226
227
228
            std::vector<FilterSharedPtr>::iterator x;
            for (x = m_filters.begin(); x != m_filters.end(); ++x)
            {
                (*x)->Initialise(m_fields, m_time);
            }
229

Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
230
            // Ensure that there is no conflict of parameters
231
            if(m_cflSafetyFactor > 0.0)
232
            {
233
                // Check final condition
234
235
236
                ASSERTL0(m_fintime == 0.0 || m_steps == 0,
                         "Final condition not unique: "
                         "fintime > 0.0 and Nsteps > 0");
237
                
238
239
240
                // Check timestep condition
                ASSERTL0(m_timestep == 0.0, 
                         "Timestep not unique: timestep > 0.0 & CFL > 0.0");
241
242
            }

Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
243
            // Check uniqueness of checkpoint output
Dave Moxey's avatar
Dave Moxey committed
244
245
246
            ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
                     (m_checktime >  0.0 && m_checksteps == 0) || 
                     (m_checktime == 0.0 && m_checksteps >  0),
247
248
                     "Only one of IO_CheckTime and IO_CheckSteps "
                     "should be set!");
249

250
251
            Timer     timer;
            bool      doCheckTime   = false;
252
            int       step          = m_initialStep;
253
254
            NekDouble intTime       = 0.0;
            NekDouble lastCheckTime = 0.0;
255
            NekDouble cpuTime       = 0.0;
Chris Cantwell's avatar
Chris Cantwell committed
256
            NekDouble elapsed       = 0.0;
257

258
259
            while (step   < m_steps ||
                   m_time < m_fintime - NekConstants::kNekZeroTol)
260
            {
261
                if (m_cflSafetyFactor)
262
                {
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
                    m_timestep = GetTimeStep(fields);
                    
                    // Ensure that the final timestep finishes at the final
                    // time, or at a prescribed IO_CheckTime.
                    if (m_time + m_timestep > m_fintime && m_fintime > 0.0)
                    {
                        m_timestep = m_fintime - m_time;
                    }
                    else if (m_checktime && 
                             m_time + m_timestep - lastCheckTime >= m_checktime)
                    {
                        lastCheckTime += m_checktime;
                        m_timestep     = lastCheckTime - m_time;
                        doCheckTime    = true;
                    }
278
                }
279
                
Chris Cantwell's avatar
Chris Cantwell committed
280
                // Perform any solver-specific pre-integration steps
281
                timer.Start();
282
283
284
285
286
                if (v_PreIntegrate(step))
                {
                    break;
                }

287
288
                fields = m_intScheme->TimeIntegrate(
                    step, m_timestep, m_intSoln, m_ode);
289
                timer.Stop();
290

291
                m_time  += m_timestep;
Chris Cantwell's avatar
Chris Cantwell committed
292
                elapsed  = timer.TimePerTest(1);
293
                intTime += elapsed;
294
                cpuTime += elapsed;
295
		
296
297
298
299
                // Write out status information
                if (m_session->GetComm()->GetRank() == 0 && 
                    !((step+1) % m_infosteps))
                {
300
301
302
303
304
305
306
307
308
                    cout << "Steps: " << setw(8)  << left << step+1 << " "
                         << "Time: "  << setw(12) << left << m_time;

                    if (m_cflSafetyFactor)
                    {
                        cout << " Time-step: " << setw(12)
                             << left << m_timestep;
                    }

309
310
                    stringstream ss;
                    ss << cpuTime << "s";
Dave Moxey's avatar
Dave Moxey committed
311
312
                    cout << " CPU Time: " << setw(8) << left
                         << ss.str() << endl;
313
                    cpuTime = 0.0;
314
                }
315

316
                // Transform data into coefficient space
317
                for (i = 0; i < nvariables; ++i)
318
                {
319
                    m_fields[m_intVariables[i]]->SetPhys(fields[i]);
320
321
322
323
                    m_fields[m_intVariables[i]]->FwdTrans_IterPerExp(
                        fields[i],
                        m_fields[m_intVariables[i]]->UpdateCoeffs());
                    m_fields[m_intVariables[i]]->SetPhysState(false);
324
                }
325

326
327
328
329
330
331
                // Perform any solver-specific post-integration steps
                if (v_PostIntegrate(step))
                {
                    break;
                }

332
                // search for NaN and quit if found
333
                int nanFound = 0;
334
335
336
337
                for (i = 0; i < nvariables; ++i)
                {
                    if (Vmath::Nnan(fields[i].num_elements(), fields[i], 1) > 0)
                    {
338
                        nanFound = 1;
339
340
                    }
                }
341
342
343
344
                m_session->GetComm()->AllReduce(nanFound,
                            LibUtilities::ReduceMax);
                ASSERTL0 (!nanFound,
                            "NaN found during time integration.");
345

Chris Cantwell's avatar
Chris Cantwell committed
346
                // Update filters
347
348
                std::vector<FilterSharedPtr>::iterator x;
                for (x = m_filters.begin(); x != m_filters.end(); ++x)
349
                {
350
                    (*x)->Update(m_fields, m_time);
351
                }
352

353
                // Write out checkpoint files
354
                if ((m_checksteps && !((step + 1) % m_checksteps)) ||
355
                     doCheckTime)
356
                {
357
358
                    if(m_HomogeneousType == eHomogeneous1D)
                    {
359
360
                        vector<bool> transformed(nfields, false);
                        for(i = 0; i < nfields; i++)
361
                        {
362
363
364
365
366
367
368
369
                            if (m_fields[i]->GetWaveSpace())
                            {
                                m_fields[i]->SetWaveSpace(false);
                                m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
                                                      m_fields[i]->UpdatePhys());
                                m_fields[i]->SetPhysState(true);
                                transformed[i] = true;
                            }
370
                        }
371
                        Checkpoint_Output(m_nchk++);
372
                        for(i = 0; i < nfields; i++)
373
                        {
374
375
376
377
                            if (transformed[i])
                            {
                                m_fields[i]->SetWaveSpace(true);
                                m_fields[i]->HomogeneousFwdTrans(
378
379
                                    m_fields[i]->GetPhys(),
                                    m_fields[i]->UpdatePhys());
380
381
                                m_fields[i]->SetPhysState(false);
                            }
382
383
384
385
                        }
                    }
                    else
                    {
386
                        Checkpoint_Output(m_nchk++);
387
                    }
388
                    doCheckTime = false;
389
                }
390

391
392
                // Step advance
                ++step;
393
            }
394
            
Chris Cantwell's avatar
Chris Cantwell committed
395
            // Print out summary statistics
396
            if (m_session->GetComm()->GetRank() == 0)
397
            {
398
399
400
401
402
                if (m_cflSafetyFactor > 0.0)
                {
                    cout << "CFL safety factor : " << m_cflSafetyFactor << endl
                         << "CFL time-step     : " << m_timestep        << endl;
                }
Chris Cantwell's avatar
Chris Cantwell committed
403

404
405
406
407
                if (m_session->GetSolverInfo("Driver") != "SteadyState")
                {
                    cout << "Time-integration  : " << intTime  << "s"   << endl;
                }
408
409
            }
            
410
            // If homogeneous, transform back into physical space if necessary.
411
412
            if(m_HomogeneousType == eHomogeneous1D)
            {
413
                for(i = 0; i < nfields; i++)
414
                {
415
416
417
418
419
420
421
                    if (m_fields[i]->GetWaveSpace())
                    {
                        m_fields[i]->SetWaveSpace(false);
                        m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
                                              m_fields[i]->UpdatePhys());
                        m_fields[i]->SetPhysState(true);
                    }
422
423
424
425
426
427
428
429
430
431
432
                }
            }
            else
            {
                for(i = 0; i < nvariables; ++i)
                {
                    m_fields[m_intVariables[i]]->SetPhys(fields[i]);
                    m_fields[m_intVariables[i]]->SetPhysState(true);
                }
            }

Chris Cantwell's avatar
Chris Cantwell committed
433
            // Finalise filters
434
435
436
            for (x = m_filters.begin(); x != m_filters.end(); ++x)
            {
                (*x)->Finalise(m_fields, m_time);
437
            }
438
            
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
439
            // Print for 1D problems
440
441
442
443
            if(m_spacedim == 1)
            {
                v_AppendOutput1D(fields);   
            }
444
        }
445
        
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
446
        /**
447
         * @brief Sets the initial conditions.
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
448
         */
449
450
        void UnsteadySystem::v_DoInitialise()
        {
451
            CheckForRestartTime(m_time);
452
            SetBoundaryConditions(m_time);
453
            SetInitialConditions(m_time);
454
        }
455
        
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
456
        /**
457
458
         * @brief Prints a summary with some information regards the 
         * time-stepping.
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
459
         */
460
        void UnsteadySystem::v_GenerateSummary(SummaryList& s)
461
        {
462
            EquationSystem::v_GenerateSummary(s);
Chris Cantwell's avatar
Chris Cantwell committed
463
464
465
466
467
468
469
470
471
472
473
            AddSummaryItem(s, "Advection",
                           m_explicitAdvection ? "explicit" : "implicit");

            if(m_session->DefinesSolverInfo("AdvectionType"))
            {
                AddSummaryItem(s, "AdvectionType",
                               m_session->GetSolverInfo("AdvectionType"));
            }

            AddSummaryItem(s, "Diffusion",
                           m_explicitDiffusion ? "explicit" : "implicit");
474

Chris Cantwell's avatar
Chris Cantwell committed
475
            if (m_session->GetSolverInfo("EQTYPE")
476
                    == "SteadyAdvectionDiffusionReaction")
477
            {
478
479
                AddSummaryItem(s, "Reaction",
                               m_explicitReaction  ? "explicit" : "implicit");
480
            }
481

482
483
484
            AddSummaryItem(s, "Time Step", m_timestep);
            AddSummaryItem(s, "No. of Steps", m_steps);
            AddSummaryItem(s, "Checkpoints (steps)", m_checksteps);
485
486
487
            AddSummaryItem(s, "Integration Type",
                           LibUtilities::TimeIntegrationMethodMap[
                               m_intScheme->GetIntegrationMethod()]);
488
489
        }
        
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
490
491
        /**
         * Stores the solution in a file for 1D problems only. This method has 
492
         * been implemented to facilitate the post-processing for 1D problems.
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
493
494
495
         */
        void UnsteadySystem::v_AppendOutput1D(
            Array<OneD, Array<OneD, NekDouble> > &solution1D)
496
        {
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
497
            // Coordinates of the quadrature points in the real physical space
498
499
500
501
502
            Array<OneD,NekDouble> x(GetNpoints());
            Array<OneD,NekDouble> y(GetNpoints());
            Array<OneD,NekDouble> z(GetNpoints());
            m_fields[0]->GetCoords(x, y, z);
            
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
503
            // Print out the solution in a txt file
504
505
506
507
            ofstream outfile;
            outfile.open("solution1D.txt");
            for(int i = 0; i < GetNpoints(); i++)
            {
508
509
                outfile << scientific << setw (17) << setprecision(16) << x[i]
                        << "  " << solution1D[0][i] << endl;
510
511
512
            }
            outfile << endl << endl;
            outfile.close();
513
514
515
516
517
518
        }

        void UnsteadySystem::v_NumericalFlux(
            Array<OneD, Array<OneD, NekDouble> > &physfield,
            Array<OneD, Array<OneD, NekDouble> > &numflux)
        {
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
519
520
            ASSERTL0(false, 
                     "This function is not implemented for this equation.");
521
522
523
524
525
526
527
        }

        void UnsteadySystem::v_NumericalFlux(
            Array<OneD, Array<OneD, NekDouble> > &physfield,
            Array<OneD, Array<OneD, NekDouble> > &numfluxX,
            Array<OneD, Array<OneD, NekDouble> > &numfluxY )
        {
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
528
529
            ASSERTL0(false, 
                     "This function is not implemented for this equation.");
530
531
532
        }

        void UnsteadySystem::v_NumFluxforScalar(
533
534
            const Array<OneD, Array<OneD, NekDouble> >               &ufield,
                  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &uflux)
535
        {
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
536
            int i, j;
537
            int nTraceNumPoints = GetTraceNpoints();
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
538
539
            int nvariables      = m_fields.num_elements();
            int nqvar           = uflux.num_elements();
540

541
542
543
544
            Array<OneD, NekDouble > Fwd     (nTraceNumPoints);
            Array<OneD, NekDouble > Bwd     (nTraceNumPoints);
            Array<OneD, NekDouble > Vn      (nTraceNumPoints, 0.0);
            Array<OneD, NekDouble > fluxtemp(nTraceNumPoints, 0.0);
545
546
547

            // Get the sign of (v \cdot n), v = an arbitrary vector

548
549
            // Evaulate upwind flux:
            // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
550
551
            for (j = 0; j < nqvar; ++j)
            {
552
                for (i = 0; i < nvariables ; ++i)
553
                {
554
555
                    // Compute Fwd and Bwd value of ufield of i direction
                    m_fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
556
557

                    // if Vn >= 0, flux = uFwd, i.e.,
558
559
                    // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
                    // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
560
561

                    // else if Vn < 0, flux = uBwd, i.e.,
562
563
                    // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
                    // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
564

565
566
                    m_fields[i]->GetTrace()->Upwind(m_traceNormals[j], 
                                                    Fwd, Bwd, fluxtemp);
567
568
569

                    // Imposing weak boundary condition with flux
                    // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
570
571
                    // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
                    // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
572
573

                    // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
574
575
                    // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
                    // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
576
577
578

                    if(m_fields[0]->GetBndCondExpansions().num_elements())
                    {
579
                        WeakPenaltyforScalar(i, ufield[i], fluxtemp);
580
581
                    }

582
583
584
585
                    // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}), 
                    // i.e,
                    // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
                    // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
586

587
588
589
590
                    // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}), 
                    // i.e,
                    // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
                    // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
591

592
593
594
595
                    Vmath::Vmul(nTraceNumPoints, 
                                m_traceNormals[j], 1, 
                                fluxtemp, 1, 
                                uflux[j][i], 1);
596
597
598
599
                }
            }
        }

600
601
        
        
602
        void UnsteadySystem::v_NumFluxforVector(
603
604
605
            const Array<OneD, Array<OneD, NekDouble> >               &ufield,
                  Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &qfield,
                  Array<OneD, Array<OneD, NekDouble> >               &qflux)
606
607
608
609
610
611
612
613
        {
            int nTraceNumPoints = GetTraceNpoints();
            int nvariables = m_fields.num_elements();
            int nqvar = qfield.num_elements();

            NekDouble C11 = 1.0;
            Array<OneD, NekDouble > Fwd(nTraceNumPoints);
            Array<OneD, NekDouble > Bwd(nTraceNumPoints);
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
614
            Array<OneD, NekDouble > Vn (nTraceNumPoints, 0.0);
615

616
617
            Array<OneD, NekDouble > qFwd     (nTraceNumPoints);
            Array<OneD, NekDouble > qBwd     (nTraceNumPoints);
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
618
            Array<OneD, NekDouble > qfluxtemp(nTraceNumPoints, 0.0);
619
620
621

            Array<OneD, NekDouble > uterm(nTraceNumPoints);

622
623
624
            // Evaulate upwind flux:
            // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
            for (int i = 0; i < nvariables; ++i)
625
            {
626
627
                qflux[i] = Array<OneD, NekDouble> (nTraceNumPoints, 0.0);
                for (int j = 0; j < nqvar; ++j)
628
                {
629
                    //  Compute Fwd and Bwd value of ufield of jth direction
630
631
632
                    m_fields[i]->GetFwdBwdTracePhys(qfield[j][i],qFwd,qBwd);

                    // if Vn >= 0, flux = uFwd, i.e.,
633
634
635
636
                    // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick 
                    // qflux = qBwd = q+
                    // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick 
                    // qflux = qBwd = q-
637
638

                    // else if Vn < 0, flux = uBwd, i.e.,
639
640
641
642
643
644
645
646
647
648
649
650
651
                    // edge::eForward, if V*n<0 <=> V*n_F<0, pick 
                    // qflux = qFwd = q-
                    // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick 
                    // qflux = qFwd = q+

                    m_fields[i]->GetTrace()->Upwind(m_traceNormals[j], 
                                                    qBwd, qFwd, 
                                                    qfluxtemp);
                    
                    Vmath::Vmul(nTraceNumPoints, 
                                m_traceNormals[j], 1, 
                                qfluxtemp, 1, 
                                qfluxtemp, 1);
652
653

                    // Generate Stability term = - C11 ( u- - u+ )
654
655
656
657
658
659
660
661
662
                    m_fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
                    
                    Vmath::Vsub(nTraceNumPoints, 
                                Fwd, 1, Bwd, 1, 
                                uterm, 1);
                    
                    Vmath::Smul(nTraceNumPoints, 
                                -1.0 * C11, uterm, 1, 
                                uterm, 1);
663

664
665
666
667
668
                    // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
                    Vmath::Vadd(nTraceNumPoints, 
                                uterm, 1, 
                                qfluxtemp, 1, 
                                qfluxtemp, 1);
669
670

                    // Imposing weak boundary condition with flux
671
                    if (m_fields[0]->GetBndCondExpansions().num_elements())
672
                    {
673
674
675
676
                        WeakPenaltyforVector(i, j, 
                                             qfield[j][i], 
                                             qfluxtemp, 
                                             C11);
677
678
679
                    }

                    // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
680
681
682
683
684
685
                    // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
                    // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
                    Vmath::Vadd(nTraceNumPoints, 
                                qfluxtemp, 1, 
                                qflux[i], 1, 
                                qflux[i], 1);
686
687
688
689
                }
            }
        }

690
691
692
693
        void UnsteadySystem::CheckForRestartTime(NekDouble &time)
        {
            if (m_session->DefinesFunction("InitialConditions"))
            {
694
                for (int i = 0; i < m_fields.num_elements(); ++i)
695
696
                {
                    LibUtilities::FunctionType vType;
697
698
699
700

                    vType = m_session->GetFunctionType(
                        "InitialConditions", m_session->GetVariable(i));

701
702
703
                    if (vType == LibUtilities::eFunctionTypeFile)
                    {
                        std::string filename
704
705
706
                            = m_session->GetFunctionFilename(
                                "InitialConditions", m_session->GetVariable(i));

707
708
709
710
711
                        fs::path pfilename(filename);

                        // redefine path for parallel file which is in directory
                        if(fs::is_directory(pfilename))
                        {
712
                            fs::path metafile("Info.xml");
713
714
715
                            fs::path fullpath = pfilename / metafile;
                            filename = LibUtilities::PortablePath(fullpath);
                        }
716
                        m_fld->ImportFieldMetaData(filename, m_fieldMetaDataMap);
717
718

                        // check to see if time defined
719
720
                        if (m_fieldMetaDataMap !=
                                LibUtilities::NullFieldMetaDataMap)
721
                        {
722
                            LibUtilities::FieldMetaDataMap::iterator iter; 
723
724
                            
                            iter = m_fieldMetaDataMap.find("Time");
725
                            if (iter != m_fieldMetaDataMap.end())
726
                            {
727
728
                                time = boost::lexical_cast<NekDouble>(
                                    iter->second);
729
730
731
732
733
734
735
736
737
                            }
                        }
                        
                        break;
                    }
                }
            }
        }
        
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
738
739
740
741
742
        void UnsteadySystem::WeakPenaltyforScalar(
            const int var,
            const Array<OneD, const NekDouble> &physfield,
                  Array<OneD,       NekDouble> &penaltyflux,
            NekDouble time)
743
        {
Dave Moxey's avatar
Dave Moxey committed
744
            int i, e, npoints, id1, id2;
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
745
            
746
747
748
749
750
751
752
753
            // Number of boundary regions
            int nbnd = m_fields[var]->GetBndCondExpansions().num_elements();
            int Nfps, numBDEdge;
            int nTraceNumPoints = GetTraceNpoints();
            int cnt = 0;

            Array<OneD, NekDouble > uplus(nTraceNumPoints);

754
755
            m_fields[var]->ExtractTracePhys(physfield, uplus);
            for (i = 0; i < nbnd; ++i)
756
757
            {
                // Number of boundary expansion related to that region
758
759
760
                numBDEdge = m_fields[var]->
                    GetBndCondExpansions()[i]->GetExpSize();
                
761
                // Evaluate boundary values g_D or g_N from input files
762
763
764
765
766
767
                LibUtilities::EquationSharedPtr ifunc = 
                    m_session->GetFunction("InitialConditions", 0);
                
                npoints = m_fields[var]->
                    GetBndCondExpansions()[i]->GetNpoints();
                
768
769
770
771
772
773
774
775
776
777
778
779
                Array<OneD,NekDouble> BDphysics(npoints);
                Array<OneD,NekDouble> x0(npoints,0.0);
                Array<OneD,NekDouble> x1(npoints,0.0);
                Array<OneD,NekDouble> x2(npoints,0.0);

                m_fields[var]->GetBndCondExpansions()[i]->GetCoords(x0,x1,x2);
                ifunc->Evaluate(x0,x1,x2,time,BDphysics);

                // Weakly impose boundary conditions by modifying flux values
                for (e = 0; e < numBDEdge ; ++e)
                {
                    // Number of points on the expansion
780
781
782
783
784
785
786
787
788
                    Nfps = m_fields[var]->
                        GetBndCondExpansions()[i]->GetExp(e)->GetNumPoints(0);
                    
                    id1 = m_fields[var]->
                        GetBndCondExpansions()[i]->GetPhys_Offset(e);
                    
                    id2 = m_fields[0]->GetTrace()->
                        GetPhys_Offset(m_fields[0]->GetTraceMap()->
                                        GetBndCondTraceToGlobalTraceMap(cnt++));
789
790

                    // For Dirichlet boundary condition: uflux = g_D
791
792
                    if (m_fields[var]->GetBndConditions()[i]->
                    GetBoundaryConditionType() == SpatialDomains::eDirichlet)
793
                    {
794
795
796
                        Vmath::Vcopy(Nfps, 
                                     &BDphysics[id1], 1, 
                                     &penaltyflux[id2], 1);
797
798
                    }
                    // For Neumann boundary condition: uflux = u+
799
800
                    else if ((m_fields[var]->GetBndConditions()[i])->
                    GetBoundaryConditionType() == SpatialDomains::eNeumann)
801
                    {
802
803
804
                        Vmath::Vcopy(Nfps, 
                                     &uplus[id2], 1, 
                                     &penaltyflux[id2], 1);
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
                    }
                }
            }
        }

        /**
         * Diffusion: Imposing weak boundary condition for q with flux
         *  uflux = g_D  on Dirichlet boundary condition
         *  uflux = u_Fwd  on Neumann boundary condition
         */
        void UnsteadySystem::WeakPenaltyforVector(
            const int var,
            const int dir,
            const Array<OneD, const NekDouble> &physfield,
            Array<OneD, NekDouble> &penaltyflux,
            NekDouble C11,
            NekDouble time)
        {
Dave Moxey's avatar
Dave Moxey committed
823
            int i, e, npoints, id1, id2;
824
825
826
827
828
829
830
831
832
            int nbnd = m_fields[var]->GetBndCondExpansions().num_elements();
            int numBDEdge, Nfps;
            int nTraceNumPoints = GetTraceNpoints();
            Array<OneD, NekDouble > uterm(nTraceNumPoints);
            Array<OneD, NekDouble > qtemp(nTraceNumPoints);
            int cnt = 0;

            m_fields[var]->ExtractTracePhys(physfield,qtemp);

833
            for (i = 0; i < nbnd; ++i)
834
            {
835
836
837
                numBDEdge = m_fields[var]->
                    GetBndCondExpansions()[i]->GetExpSize();
                
838
                // Evaluate boundary values g_D or g_N from input files
839
840
841
842
843
                LibUtilities::EquationSharedPtr ifunc = 
                    m_session->GetFunction("InitialConditions", 0);
                
                npoints = m_fields[var]->
                    GetBndCondExpansions()[i]->GetNpoints();
844
845
846
847
848
849
850
851
852
853
854
855

                Array<OneD,NekDouble> BDphysics(npoints);
                Array<OneD,NekDouble> x0(npoints,0.0);
                Array<OneD,NekDouble> x1(npoints,0.0);
                Array<OneD,NekDouble> x2(npoints,0.0);

                m_fields[var]->GetBndCondExpansions()[i]->GetCoords(x0,x1,x2);
                ifunc->Evaluate(x0,x1,x2,time,BDphysics);

                // Weakly impose boundary conditions by modifying flux values
                for (e = 0; e < numBDEdge ; ++e)
                {
856
857
                    Nfps = m_fields[var]->
                        GetBndCondExpansions()[i]->GetExp(e)->GetNumPoints(0);
858

859
860
861
862
863
864
865
866
867
868
869
                    id1 = m_fields[var]->
                        GetBndCondExpansions()[i]->GetPhys_Offset(e);
                    
                    id2 = m_fields[0]->GetTrace()->
                        GetPhys_Offset(m_fields[0]->GetTraceMap()->
                                       GetBndCondTraceToGlobalTraceMap(cnt++));

                    // For Dirichlet boundary condition: 
                    //qflux = q+ - C_11 (u+ -    g_D) (nx, ny)
                    if(m_fields[var]->GetBndConditions()[i]->
                    GetBoundaryConditionType() == SpatialDomains::eDirichlet)
870
                    {
871
872
873
874
                        Vmath::Vmul(Nfps, 
                                    &m_traceNormals[dir][id2], 1, 
                                    &qtemp[id2], 1, 
                                    &penaltyflux[id2], 1);
875
876
                    }
                    // For Neumann boundary condition: qflux = g_N
877
878
                    else if((m_fields[var]->GetBndConditions()[i])->
                    GetBoundaryConditionType() == SpatialDomains::eNeumann)
879
                    {
880
881
882
883
                        Vmath::Vmul(Nfps,
                                    &m_traceNormals[dir][id2], 1, 
                                    &BDphysics[id1], 1, 
                                    &penaltyflux[id2], 1);
884
885
886
887
888
                    }
                }
            }
        }
	
889
        /**
890
891
892
893
894
895
         * @brief Return the timestep to be used for the next step in the
         * time-marching loop.
         *
         * This function can be overloaded to facilitate solver which utilise a
         * CFL (or other) parameter to determine a maximum timestep under which
         * the problem remains stable.
896
         */
897
898
        NekDouble UnsteadySystem::GetTimeStep(
            const Array<OneD, const Array<OneD, NekDouble> > &inarray)
899
        {
900
            return v_GetTimeStep(inarray);
901
902
903
        }
        
        /**
904
905
906
907
         * @brief Return the timestep to be used for the next step in the
         * time-marching loop.
         *
         * @see UnsteadySystem::GetTimeStep
908
         */
909
910
        NekDouble UnsteadySystem::v_GetTimeStep(
            const Array<OneD, const Array<OneD, NekDouble> > &inarray)
911
912
        {
            ASSERTL0(false, "Not defined for this class");
913
            return 0.0;
914
        }
915
916
917
918
919
920
921
922
923
924

        bool UnsteadySystem::v_PreIntegrate(int step)
        {
            return false;
        }

        bool UnsteadySystem::v_PostIntegrate(int step)
        {
            return false;
        }
Dave Moxey's avatar
Dave Moxey committed
925

926
927
928
929
        bool UnsteadySystem::v_SteadyStateCheck(int step)
        {
            return false;
        }
930

Chris Cantwell's avatar
Chris Cantwell committed
931
932
933
        void UnsteadySystem::SVVVarDiffCoeff(
            const Array<OneD, Array<OneD, NekDouble> >  vel,
                  StdRegions::VarCoeffMap              &varCoeffMap)
934
935
936
937
938
        {
            int phystot = m_fields[0]->GetTotPoints();
            int nvel = vel.num_elements();

            Array<OneD, NekDouble> varcoeff(phystot),tmp;
Chris Cantwell's avatar
Chris Cantwell committed
939

940
941
942
943
944
945
946
947
948
949
950
951
            // calculate magnitude of v
            Vmath::Vmul(phystot,vel[0],1,vel[0],1,varcoeff,1);
            for(int n = 1; n < nvel; ++n)
            {
                Vmath::Vvtvp(phystot,vel[n],1,vel[n],1,varcoeff,1,varcoeff,1);
            }
            Vmath::Vsqrt(phystot,varcoeff,1,varcoeff,1);

            for(int i = 0; i < m_fields[0]->GetNumElmts(); ++i)
            {
                int offset = m_fields[0]->GetPhys_Offset(i);
                int nq = m_fields[0]->GetExp(i)->GetTotPoints();
Chris Cantwell's avatar
Chris Cantwell committed
952
953
954
                Array<OneD, NekDouble> unit(nq,1.0);

                int nmodes = 0;
955
956
957

                for(int n = 0; n < m_fields[0]->GetExp(i)->GetNumBases(); ++n)
                {
Chris Cantwell's avatar
Chris Cantwell committed
958
959
                    nmodes = max(nmodes,
                                 m_fields[0]->GetExp(i)->GetBasisNumModes(n));
960
                }
Chris Cantwell's avatar
Chris Cantwell committed
961

962
963
964
965
966
967
                NekDouble h = m_fields[0]->GetExp(i)->Integral(unit);
                h = pow(h,(NekDouble) (1.0/nvel))/((NekDouble) nmodes);

                Vmath::Smul(nq,h,varcoeff+offset,1,tmp = varcoeff+offset,1);
            }

Chris Cantwell's avatar
Chris Cantwell committed
968
969
            // set up map with eVarCoffLaplacian key
            varCoeffMap[StdRegions::eVarCoeffLaplacian] = varcoeff;
970
        }
971
972
    }
}