TimeIntegrationDemo.cpp 23 KB
Newer Older
Peter Vos's avatar
Peter Vos committed
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
39
40
41
///////////////////////////////////////////////////////////////////////////////
//
// File TimeIntegrationDemo.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Scientific Computing and Imaging Institute,
// University of Utah (USA) and Department of Aeronautics, Imperial
// College London (UK).
//
// 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:
//
///////////////////////////////////////////////////////////////////////////////

//--------------------------------------------------
// This file illustrates the use of the time-stepping framework.
//
// For more information, please consult the following reference:
//
// Vos, P.E.J., Eskilsson, C., Bolis, A., Chun, S., Kirby, R.M. and Sherwin, S.J.
42
// "A Generic Framework for Time-Stepping PDEs: general linear methods,
Peter Vos's avatar
Peter Vos committed
43
44
45
46
47
48
//  object-oriented implementation and application to fluid problems"
// International Journal of Computational Fluid Dynamics, to appear

// It solves the one-dimensional advection-diffusion problem, defined as
//
//  |    du     du     d^2 u
49
//  |    -- + V -- = D -----,
Peter Vos's avatar
Peter Vos committed
50
51
52
53
54
55
//  |    dt     dx     d x^2
//  |
//  |    subject to:
//  |    - periodic boundary conditions
//  |    - the initial condition
//  |        u(x,0) = sin(2*pi*k*x)
56
57
//  |
//  |    and with  x = [0,1]
Peter Vos's avatar
Peter Vos committed
58
59
60
61
//  |              t = [0,1]
//  |              U = 1
//  |              D = 0.05
//  |              k = 1
62
//  |
Peter Vos's avatar
Peter Vos committed
63
64
65
66
67
68
//
// using the finite difference method.
// The exact solution of the equation above is given by
//
//  u(x,t) = exp(-D * (2*pi*k)^2 * t) * sin(2*pi*k * (x - U*t) )
//
69
// The output is written out to the files
Peter Vos's avatar
Peter Vos committed
70
71
72
73
74
75
76
77
78
79
//
//   - OneDfinDiffAdvDiffSolverOutput.dat (containing the data)
//   - OneDfinDiffAdvDiffSolverOutput.p   (containing a gnuplot script)
//
// and can be visualised by gnuplot using the command
//
//    gnuplot OneDfinDiffAdvDiffSolverOutput.p
//
//--------------------------------------------------
#include <fstream>
80
#include <iostream>
Peter Vos's avatar
Peter Vos committed
81
#include <iomanip>
82
83
84

#include <boost/program_options.hpp>
#include <boost/algorithm/string.hpp>
Peter Vos's avatar
Peter Vos committed
85
86
87
88
89
90
91
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/TimeIntegration/TimeIntegrationScheme.h>

using namespace std;
using namespace Nektar;
using namespace Nektar::LibUtilities;

92
93
namespace po = boost::program_options;

94
// We first implement a class that represents
Peter Vos's avatar
Peter Vos committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
// the 1D finite difference solver
class OneDfinDiffAdvDiffSolver
{
public:
    // constructor based upon the discretisation details
    OneDfinDiffAdvDiffSolver(int nPoints, int nTimeSteps):
        m_x0(0.0),
        m_xend(1.0),
        m_nPoints(nPoints),
        m_dx((m_xend-m_x0)/((double) m_nPoints-1.0)),
        m_t0(0.0),
        m_tend(1.0),
        m_nTimeSteps(nTimeSteps),
        m_dt((m_tend-m_t0)/(double) m_nTimeSteps),
        m_wavenumber(1.0),
        m_U(1.0),
        m_D(0.05)
    {
    }

    // -----------------------------------------------------------------
    // ---- These functions/methods below are the routines which will be used by
    // ---- the TimeIntegration framework (and are required for using it)...
    // ---- The implementation of these functions can be found at the end of the file
    void HelmSolve(const Array<OneD, const Array<OneD, double> >& inarray,
                         Array<OneD,       Array<OneD, double> >& outarray,
                   const NekDouble time,
122
123
                   const NekDouble lambda) const;

Peter Vos's avatar
Peter Vos committed
124
125
126
    void EvaluateAdvectionTerm(const Array<OneD, const  Array<OneD, double> >& inarray,
                                     Array<OneD,        Array<OneD, double> >& outarray,
                               const NekDouble time) const;
127

Peter Vos's avatar
Peter Vos committed
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
    void Project(const Array<OneD, const  Array<OneD, double> >& inarray,
                       Array<OneD,        Array<OneD, double> >& outarray,
                 const NekDouble time) const;
    // -----------------------------------------------------------------


    // -----------------------------------------------------------------
    // Below are some additional functions
    int GetNpoints() const;

    void EvaluateExactSolution(Array<OneD, Array<OneD, double> >& outarray,
                               const NekDouble time) const;

    double EvaluateL2Error(const Array<OneD, const  Array<OneD, double> >& approx,
                           const Array<OneD, const  Array<OneD, double> >& exact) const;

    void AppendOutput(ofstream& outfile,
                      const Array<OneD, const  Array<OneD, double> >& approx,
                      const Array<OneD, const  Array<OneD, double> >& exact) const;

    void GenerateGnuplotScript() const;

    double GetInitialTime() const;

    double GetTimeStep() const;
    // -----------------------------------------------------------------

private:
    // spatial discretisation
    double m_x0;         // the left boundary of the domain
    double m_xend;       // the right boundary of the domain
    int    m_nPoints;    // the number of grid-points used in the finite difference method
    double m_dx;         // the distance between 2 grid points

    // temporal discretisation
    double m_t0;         // the initial time
    double m_tend;       // the end time
    int    m_nTimeSteps; // the number of time-steps
    double m_dt;         // the size of a time-step

    // value of the coefficients
169
    double m_wavenumber; // wave number
Peter Vos's avatar
Peter Vos committed
170
171
172
    double m_U;          // advection speed
    double m_D;          // diffusion coefficient

173
174

    void solveTriDiagMatrix (int n, double alpha, double beta,
Peter Vos's avatar
Peter Vos committed
175
                             const Array<OneD, const double>& inarray,
176
                             Array<OneD,       double>& outarray) const;
Peter Vos's avatar
Peter Vos committed
177
178
};

179

Peter Vos's avatar
Peter Vos committed
180
int main(int argc, char *argv[])
181
{
182
183
184
185
186
187
188
189
190
191
192
    po::options_description desc("Available options");
    desc.add_options()
        ("help, h",         "Produce this help message.")
        ("Npoints, np",     po::value<int>(),
                            "the number of grid points to be used.")
        ("Ntimesteps, nt",  po::value<int>(),
                            "the number of timesteps to be used.")
        ("NTimeIntegrationMethod, nm",  po::value<int>(),
                            "TimeIntegrationMethod is a number in the range [1,5]."
                            "and defines the time-integration method to be used, i.e"
                            "- 1: 1st order multi-step IMEX scheme (Euler Backwards/Euler Forwards)"
193
194
195
                            "- 2: 2nd order multi-step IMEX scheme"
                            "- 3: 3rd order multi-step IMEX scheme"
                            "- 4: 2nd order multi-stage DIRK IMEX scheme"
196
                            "- 5: 3nd order multi-stage DIRK IMEX scheme"
197
                            "- 6: 2nd order IMEX Gear (Extrapolated Gear/SBDF-2)");
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
    po::variables_map vm;
    try
    {
        po::store(po::command_line_parser(argc, argv).
                  options(desc).run(), vm);
        po::notify(vm);
    }
    catch (const std::exception& e)
    {
        cerr << e.what() << endl;
        cerr << desc;
        return 1;
    }

    if (!vm.count("Npoints") || !vm.count("Ntimesteps") || !vm.count("NTimeIntegrationMethod"))
    {
214
        cerr << "Usage: Project1D --Npoints nPoints --Ntimesteps nTimesteps --TimeIntegrationMethod nMethod" << endl;
215
216
        cerr << "Where  - Npoints is the number of grid points to be used" << endl;
        cerr << "         for the finite difference discretisation" << endl;
217
        cerr << "       - Ntimesteps is the number of timesteps to be used" << endl;
218
219
220
221
222
223
224
225
        cerr << "         for the time-integration method" << endl;
        cerr << "       - TimeIntegrationMethod is a number in the range [1,8]" << endl;
        cerr << "         and defines the time-integration method to be used, i.e." << endl;
        cerr << "           - 1: 1st order multi-step IMEX scheme (Euler Backwards/Euler Forwards)" << endl;
        cerr << "           - 2: 2nd order multi-step IMEX scheme" << endl;
        cerr << "           - 3: 3rd order multi-step IMEX scheme" << endl;
        cerr << "           - 4: 2nd order multi-stage DIRK IMEX scheme" << endl;
        cerr << "           - 5: 3nd order multi-stage DIRK IMEX scheme" << endl;
226
227
228
229
230
231
232
233
234
235
        cerr << "           - 6: 2nd order IMEX Gear (Extrapolated Gear/SBDF-2)" << endl;
        cerr << "           - 7: 2nd order Crank-Nicolson/Adams-Bashforth (CNAB)" << endl;
        cerr << "           - 8: 2nd order Modified Crank-Nicolson/Adams-Bashforth (MCNAB)" << endl;
        return 1;
    }

    int nPoints = vm["Npoints"].as<int>();
    int nTimesteps = vm["Ntimesteps"].as<int>();
    int nMethod = vm["NTimeIntegrationMethod"].as<int>();

Peter Vos's avatar
Peter Vos committed
236
237
238
239
    // Open a file for writing the solution
    ofstream outfile;
    outfile.open("OneDfinDiffAdvDiffSolverOutput.dat");

240

Peter Vos's avatar
Peter Vos committed
241
242
243
244
245
246
247

    // -----------------------------------------------------------------------------
    // THE IMPLEMENTATION BELOW SHOWS HOW THE TIME-STEPPING FRAMEWORK CAN BE
    // USED FOR TIME-INTEGRATION PDEs

    // 1. THE SPATIAL DISCRETISATION
    //    Create an object of the OneDfinDiffAdvDiffSolver class.
248
    //    This class can be thought of as representing the
Peter Vos's avatar
Peter Vos committed
249
250
    //    spatial (finite difference) discretisation.
    OneDfinDiffAdvDiffSolver* solver = new OneDfinDiffAdvDiffSolver(nPoints,nTimesteps);
251
    //    After this spatial discretisation, the PDE has actually been
Peter Vos's avatar
Peter Vos committed
252
253
254
    //    reduced (through the method-of-lines) to an ODE.
    //    In order to use the time-stepping framework, we need to give it the necessary
    //    information about this ODE.
255
    //    Therefore, we create an oject of the class TimeIntegrationSchemeOperators that
Peter Vos's avatar
Peter Vos committed
256
257
258
259
260
261
262
263
264
265
266
267
    //    contains a 'function pointer' (in fact a 'functor') to the
    //    - explicit term of the ODE (i.e. the advection term)
    //    - implicit solve routine (i.e. the Helmholtz solver)
    //    - projection operator (i.e. the identity operator in this case)
    LibUtilities::TimeIntegrationSchemeOperators ode;
    ode.DefineOdeRhs        (&OneDfinDiffAdvDiffSolver::EvaluateAdvectionTerm,solver);
    ode.DefineImplicitSolve (&OneDfinDiffAdvDiffSolver::HelmSolve,            solver);
    ode.DefineProjection    (&OneDfinDiffAdvDiffSolver::Project,              solver);

    // 2. THE TEMPORAL DISCRETISATION
    // 2.1 Read in which method should be used.
    //     For a multi-step scheme, also set up
268
    //     which method should be used for appropriately
Peter Vos's avatar
Peter Vos committed
269
270
271
    //     starting up the system
    Array<OneD, TimeIntegrationMethod> method;
    int nSteps = 1;
Chris Cantwell's avatar
Chris Cantwell committed
272
    switch (nMethod)
Peter Vos's avatar
Peter Vos committed
273
274
275
276
277
278
    {
    case 1 :
        {
            nSteps = 1;
            method = Array<OneD, TimeIntegrationMethod>(nSteps);
            method[0] = eIMEXOrder1;
279
        }
Peter Vos's avatar
Peter Vos committed
280
281
282
283
284
285
        break;
    case 2 :
        {
            nSteps = 2;
            method = Array<OneD, TimeIntegrationMethod>(nSteps);
            method[0] = eIMEXdirk_2_3_2; // the start-up method for step 1
286
287
            method[1] = eIMEXOrder2;
        }
Peter Vos's avatar
Peter Vos committed
288
289
290
291
292
293
294
295
        break;
    case 3 :
        {
            nSteps = 3;
            method = Array<OneD, TimeIntegrationMethod>(nSteps);
            method[0] = eIMEXdirk_3_4_3; // the start-up method for step 1
            method[1] = eIMEXdirk_3_4_3; // the start-up method for step 2
            method[2] = eIMEXOrder3;
296
        }
Peter Vos's avatar
Peter Vos committed
297
298
299
300
301
302
        break;
    case 4 :
        {
            nSteps = 1;
            method = Array<OneD, TimeIntegrationMethod>(nSteps);
            method[0] = eIMEXdirk_2_3_2;
303
        }
Peter Vos's avatar
Peter Vos committed
304
305
306
307
308
309
        break;
    case 5 :
        {
            nSteps = 1;
            method = Array<OneD, TimeIntegrationMethod>(nSteps);
            method[0] = eIMEXdirk_3_4_3;
310
        }
Peter Vos's avatar
Peter Vos committed
311
        break;
312
313
314
315
316
317
318
    case 6 :
        {
            nSteps = 2;
            method = Array<OneD, TimeIntegrationMethod>(nSteps);
            method[0] = eIMEXdirk_2_2_2;
            method[1] = eIMEXGear;

319
        }
320
321
322
323
324
325
326
327
328
        break;
    case 7 :
        {
            nSteps = 3;
            method = Array<OneD, TimeIntegrationMethod>(nSteps);
            method[0] = eIMEXdirk_3_4_3;
            method[1] = eIMEXdirk_3_4_3;
            method[2] = eCNAB;

329
        }
330
331
332
333
334
335
336
337
338
        break;
    case 8 :
        {
            nSteps = 3;
            method = Array<OneD, TimeIntegrationMethod>(nSteps);
            method[0] = eIMEXdirk_3_4_3;
            method[1] = eIMEXdirk_3_4_3;
            method[2] = eMCNAB;

339
        }
340
        break;
341
342
343
344
345
346
347
348
349
350
    default :
        {
            cerr << "The third argument defines the time-integration method to be used" << endl;
            cerr << "and should be a number in the range [1,6], i.e." << endl;
            cerr << "  - 1: 1st order multi-step IMEX scheme (Euler Backwards/Euler Forwards)" << endl;
            cerr << "  - 2: 2nd order multi-step IMEX scheme" << endl;
            cerr << "  - 3: 3rd order multi-step IMEX scheme" << endl;
            cerr << "  - 4: 2nd order multi-stage DIRK IMEX scheme" << endl;
            cerr << "  - 5: 3rd order multi-stage DIRK IMEX scheme" << endl;
            cerr << "  - 6: 2nd order IMEX Gear (Extrapolated Gear/SBDF-2)" << endl;
351
352
            cerr << "  - 7: 2nd order Crank-Nicolson/Adams-Bashforth (CNAB)" << endl;
            cerr << "  - 8: 2nd order Modified Crank-Nicolson/Adams-Bashforth (MCNAB)" << endl;
353
            exit(1);
Peter Vos's avatar
Peter Vos committed
354
355
356
357
358
359
360
361
362
363
        }
    }


    // 2.2 Create objects of the time-integration framework.
    //     These can later be used for actually doing the time-integration
    Array<OneD, LibUtilities::TimeIntegrationSchemeSharedPtr> IntScheme(nSteps);
    for(int i = 0; i < nSteps; i++)
    {
        TimeIntegrationSchemeKey IntKey(method[i]);
364
        IntScheme[i] = LibUtilities::TimeIntegrationSchemeManager()[IntKey];
Peter Vos's avatar
Peter Vos committed
365
366
    }

367
    // 2.3 Initialise some arrays that contain the numerical and
Peter Vos's avatar
Peter Vos committed
368
369
370
371
372
373
374
375
376
377
378
379
    //     analytical solutions
    double t0 = solver->GetInitialTime();
    Array<OneD, Array<OneD, double> > fidifsol(1); // Array containing the numerical solution
    Array<OneD, Array<OneD, double> > exactsol(1); // Array containing the exact solution
    fidifsol[0] = Array<OneD, double>(nPoints);
    exactsol[0] = Array<OneD, double>(nPoints);
    solver->EvaluateExactSolution(fidifsol,t0);  // Set the initial condition
    solver->EvaluateExactSolution(exactsol,t0);  // Set the initial condition
    solver->AppendOutput(outfile,fidifsol,exactsol); // Write the initial condition to a file

    // 2.4 Initialize the time-integration scheme
    double dt = solver->GetTimeStep();
380
381
382
    LibUtilities::TimeIntegrationSolutionSharedPtr sol;
    sol = IntScheme[nSteps-1]->InitializeScheme(dt,fidifsol,t0,ode);

Peter Vos's avatar
Peter Vos committed
383
    // 2.5 Do the time-integration
384
    double t  = t0;
Peter Vos's avatar
Peter Vos committed
385
386
    int whichscheme;
    for(int i = 0; i < nTimesteps; i++)
387
    {
Peter Vos's avatar
Peter Vos committed
388
389
        t = t0 + i*dt;

390
        whichscheme = (i<nSteps)?i:(nSteps-1); // For multi-step schemes,
Peter Vos's avatar
Peter Vos committed
391
392
393
394
395
396
                                               // the first steps should use the start-up scheme
        fidifsol = IntScheme[whichscheme]->TimeIntegrate(dt,sol,ode); // Time-integration for 1 time-step


        solver->EvaluateExactSolution(exactsol,t);       // Calculate the exact solution
        solver->AppendOutput(outfile,fidifsol,exactsol); // Dump the output to a file
397
398
    }

Peter Vos's avatar
Peter Vos committed
399
    // Calculate the error and dump to screen
400
    cout << "L 2 error :" <<  solver->EvaluateL2Error(fidifsol,exactsol) << endl;
Peter Vos's avatar
Peter Vos committed
401
402
403
404
405

    // Some more writing out the results
    solver->GenerateGnuplotScript();
    outfile.close();

406

Peter Vos's avatar
Peter Vos committed
407
408
    delete solver;

409
  return 0;
Peter Vos's avatar
Peter Vos committed
410
411
412
413
414
415
416
417
}

void OneDfinDiffAdvDiffSolver::HelmSolve(const Array<OneD, const Array<OneD, double> >& inarray,
                                               Array<OneD,       Array<OneD, double> >& outarray,
                                         const NekDouble time,
                                         const NekDouble lambda) const
{
    // This function implements a 1D finite difference helmholtz solver.
418
    // The 1D Helmholtz equation leads to a cyclic triadiagonal matrix to be
Peter Vos's avatar
Peter Vos committed
419
420
421
422
423
424
425
426
    // solved.
    // In this function, this is solved as follows:
    // - Based upon the static condensation technique, we first solve
    //   for the periodic grid-points (i.e. the boundary nodes)
    // - Next, we solve for the interior grid-points. The associated tridiagonal
    //   system is solved based upon the Thomas algorithm
    double a = - m_D * lambda / (m_dx * m_dx); //off diagonal term
    double b  = 1.0 + 2.0 * lambda * m_D / (m_dx * m_dx); // diagonal term of triadiagonal matrix
427

Peter Vos's avatar
Peter Vos committed
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
    int nIntPoints = m_nPoints-2;

    Array<OneD, double> invD_f(nIntPoints);
    solveTriDiagMatrix(nIntPoints,a,b,inarray[0]+1,invD_f);

    Array<OneD, double> C(nIntPoints,0.0);
    Array<OneD, double> invD_C(nIntPoints,0.0);
    C[0]     = a;
    C[nIntPoints-1] = a;
    solveTriDiagMatrix(nIntPoints,a,b,C,invD_C);

    outarray[0][0] = ( inarray[0][0] - a*(invD_f[0]+invD_f[nIntPoints-1]) ) / (b - a*(invD_C[0]+invD_C[nIntPoints-1]) );
    outarray[0][m_nPoints-1] = outarray[0][0];

    Array<OneD, double> f(nIntPoints);
    for(int i = 0; i < nIntPoints; i++)
    {
        f[i] = inarray[0][i+1];
    }
    f[0] -= outarray[0][0]*a;
    f[nIntPoints-1] -= outarray[0][0]*a;

    Array<OneD, double> tmp;
    solveTriDiagMatrix(nIntPoints,a,b,f,tmp = outarray[0]+1); // Calls the Thomas algorithm
}

void OneDfinDiffAdvDiffSolver::EvaluateAdvectionTerm(const Array<OneD, const  Array<OneD, double> >& inarray,
                                                           Array<OneD,        Array<OneD, double> >& outarray,
                                                     const NekDouble time) const
{
    // The advection term can be evaluated using central or upwind differences
    if(true)
    {
        // central differences
        outarray[0][0]           = - m_U * (inarray[0][1]-inarray[0][m_nPoints-2]) / (2.0 * m_dx);
        outarray[0][m_nPoints-1] = outarray[0][0];
464

Peter Vos's avatar
Peter Vos committed
465
466
        for(int i = 1; i < m_nPoints-1; i++)
        {
467
            outarray[0][i] = - m_U * (inarray[0][i+1]-inarray[0][i-1]) / (2.0 * m_dx);
Peter Vos's avatar
Peter Vos committed
468
469
470
471
        }
    }
    else
    {
472
        // upwind differences
Peter Vos's avatar
Peter Vos committed
473
474
475
476
477
478
        for(int i = 1; i < m_nPoints; i++)
        {
            outarray[0][i] = - m_U * (inarray[0][i]-inarray[0][i-1]) / (m_dx);
        }
        outarray[0][0] = outarray[0][m_nPoints-1];

479
    }
Peter Vos's avatar
Peter Vos committed
480
481
482
483
484
485
486
487
488
489
490
491
492
493

}

void OneDfinDiffAdvDiffSolver::Project(const Array<OneD, const  Array<OneD, double> >& inarray,
                                             Array<OneD,        Array<OneD, double> >& outarray,
                                       const NekDouble time) const
{
    // This is simply the identity operator for this case
    for(int i = 0; i < m_nPoints; i++)
    {
        outarray[0][i] = inarray[0][i];
    }
}

494
void OneDfinDiffAdvDiffSolver::solveTriDiagMatrix (int n, double a, double b,
Peter Vos's avatar
Peter Vos committed
495
496
497
498
499
500
501
502
                                                   const Array<OneD, const double>& inarray,
                                                         Array<OneD,       double>& outarray) const
{
    // Implementation of the Thomas algorithm for Tridiaginol systems
    Array<OneD, double> cprime(n);
    Array<OneD, double> dprime(n);

    cprime[0] = a / b;
503
    dprime[0] = inarray[0] / b;
Peter Vos's avatar
Peter Vos committed
504
505
506
507
508
509
510
511
512

    for(int i = 1; i < n; i++)
    {
        double id = 1.0 / (b - cprime[i-1] * a);
        cprime[i] = a * id;
        dprime[i] = (inarray[i] - dprime[i-1] * a) * id;
    }

    outarray[n-1] = dprime[n-1];
513
    for (int i = n - 2; i >= 0; i--)
Peter Vos's avatar
Peter Vos committed
514
515
516
517
518
519
520
521
    {
        outarray[i] = dprime[i] - cprime[i] * outarray[i+1];
    }

}
int OneDfinDiffAdvDiffSolver::GetNpoints() const
{
    return m_nPoints;
522
}
Peter Vos's avatar
Peter Vos committed
523
524
525
526
527
528
529
530

void OneDfinDiffAdvDiffSolver::EvaluateExactSolution(Array<OneD, Array<OneD, double> >& outarray,
                                                     const NekDouble time) const
{
    double x;
    for(int i = 0; i < m_nPoints; i++)
    {
        x = m_x0 + i*m_dx;
531
        outarray[0][i] = exp(-m_D * 2.0 * 2.0 * M_PI * M_PI * m_wavenumber*m_wavenumber*time) *
Peter Vos's avatar
Peter Vos committed
532
533
534
535
536
            sin( 2.0 * m_wavenumber * M_PI * (x - m_U*time) );
    }
}
double OneDfinDiffAdvDiffSolver::EvaluateL2Error(const Array<OneD, const  Array<OneD, double> >& approx,
                                                 const Array<OneD, const  Array<OneD, double> >& exact) const
537
{
Peter Vos's avatar
Peter Vos committed
538
539
    double a = 0.0;
    double b = 0.0;
540

Peter Vos's avatar
Peter Vos committed
541
542
543
544
545
    for(int i = 0; i < m_nPoints; i++)
    {
        a += (approx[0][i]-exact[0][i])*(approx[0][i]-exact[0][i]);
        b += exact[0][i]*exact[0][i];
    }
546

Peter Vos's avatar
Peter Vos committed
547
548
549
550
551
552
553
554
555
    return sqrt(a/b);
}

void OneDfinDiffAdvDiffSolver::AppendOutput(ofstream& outfile,
                                            const Array<OneD, const  Array<OneD, double> >& approx,
                                            const Array<OneD, const  Array<OneD, double> >& exact) const
{
    for(int i = 0; i < m_nPoints; i++)
    {
556
557
558
        outfile << scientific
                << setw (17)
                << setprecision(10)
Peter Vos's avatar
Peter Vos committed
559
                << m_x0 + i*m_dx
560
561
562
563
                << "  "
                << approx[0][i]
                << "  "
                << exact[0][i]
Peter Vos's avatar
Peter Vos committed
564
565
566
567
568
569
570
571
572
                << endl;
    }
    outfile << endl << endl;
}

void OneDfinDiffAdvDiffSolver::GenerateGnuplotScript() const
{
    ofstream outfile;
    outfile.open("OneDfinDiffAdvDiffSolverOutput.p");
573

Peter Vos's avatar
Peter Vos committed
574
    outfile << "# Gnuplot script file" << endl;
575
576
577
578
579
    outfile << "set   autoscale" << endl;
    outfile << "unset log" << endl;
    outfile << "unset label" << endl;
    outfile << "set xtic auto" << endl;
    outfile << "set ytic auto" << endl;
Peter Vos's avatar
Peter Vos committed
580
581
582
583
584
    outfile << "set title \"Finite Difference Solution to the 1D advection-diffusion equation\"" << endl;
    outfile << "set xlabel \"x\"" << endl;
    outfile << "set ylabel \"u\"" << endl;
    outfile << "set xr [" << m_x0 << ":" << m_xend << "]" << endl;
    outfile << "set yr [-1.0:1.0]" << endl;
585

Peter Vos's avatar
Peter Vos committed
586
587
588
589
590
591
592
593
594
595
596
597
    double t;
    for(int i=0; i <= m_nTimeSteps; i++)
    {
        t = m_t0+i*m_dt;
        outfile << "plot    \"OneDfinDiffAdvDiffSolverOutput.dat\" ";
        outfile << "using 1:2 index ";
        outfile << i << " title 'Finite Difference Solution (t=" << t << ")' with linespoints , ";
        outfile << "\"OneDfinDiffAdvDiffSolverOutput.dat\" ";
        outfile << "using 1:3 index ";
        outfile << i << " title 'Exact Solution (t=" << t << ")' with linespoints" << endl;
        outfile << "pause " << 4.0/m_nTimeSteps << endl;
    }
598

Peter Vos's avatar
Peter Vos committed
599
600
601
602
603
604
    outfile.close();
}

double OneDfinDiffAdvDiffSolver::GetInitialTime() const
{
    return m_t0;
605
}
Peter Vos's avatar
Peter Vos committed
606
607
608
609
610
611

double OneDfinDiffAdvDiffSolver::GetTimeStep() const
{
    return m_dt;
}