ExpListHomogeneous1D.cpp 50.3 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 ExpListHomogeneous1D.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.
//
32
// Description: An ExpList which is homogeneous in 1-direction
33
34
35
//
///////////////////////////////////////////////////////////////////////////////

36
#include <MultiRegions/ExpListHomogeneous1D.h>
37
#include <LibUtilities/Foundations/ManagerAccess.h>  // for PointsManager, etc
38
#include <StdRegions/StdSegExp.h>
39
#include <StdRegions/StdPointExp.h>
40
#include <LocalRegions/Expansion.h>
41
#include <LocalRegions/Expansion2D.h>
42

43
44
using namespace std;

45
46
47
48
49
50
51
52
53
54
55
56
57
namespace Nektar
{
    namespace MultiRegions
    {
        // Forward declaration for typedefs
        ExpListHomogeneous1D::ExpListHomogeneous1D():
            ExpList(),
            m_homogeneousBasis(LibUtilities::NullBasisSharedPtr),
            m_lhom(1),
            m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr())
        {
        }

Chris Cantwell's avatar
Chris Cantwell committed
58
        ExpListHomogeneous1D::ExpListHomogeneous1D(const LibUtilities::SessionReaderSharedPtr
59
                &pSession,const LibUtilities::BasisKey &HomoBasis, const NekDouble lhom, const bool useFFT, const bool dealiasing):
60
            ExpList(pSession),
61
            m_useFFT(useFFT),
62
63
64
            m_lhom(lhom),
            m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr()),
            m_dealiasing(dealiasing)
65
        {
66
            ASSERTL2(HomoBasis != LibUtilities::NullBasisKey,"Homogeneous Basis is a null basis");
67
            
Yan Bao's avatar
Yan Bao committed
68
            m_homogeneousBasis = LibUtilities::BasisManager()[HomoBasis];
69

70
71
72
73
74
            m_StripZcomm = m_session->DefinesSolverInfo("HomoStrip") ?
                           m_comm->GetColumnComm()->GetColumnComm()  :
                           m_comm->GetColumnComm();

            m_transposition = MemoryManager<LibUtilities::Transposition>
75
                                ::AllocateSharedPtr(HomoBasis, m_comm, m_StripZcomm);
76
77
78
79

            m_planes = Array<OneD,ExpListSharedPtr>(
                                m_homogeneousBasis->GetNumPoints() /
                                m_StripZcomm->GetSize());
80

81
82
83
84
85
86
87
88
89
90
91
92
93
            ASSERTL0( m_homogeneousBasis->GetNumPoints() %
                      m_StripZcomm->GetSize() == 0,
                      "HomModesZ should be a multiple of npz.");

            if (  (m_homogeneousBasis->GetBasisType() !=
                    LibUtilities::eFourierHalfModeRe)
               && (m_homogeneousBasis->GetBasisType() !=
                    LibUtilities::eFourierHalfModeIm) )
            {
                ASSERTL0( m_planes.num_elements() % 2 == 0,
                      "HomModesZ/npz should be an even integer.");
            }

94
95
            if(m_useFFT)
            {
96
97
                m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance(
                                "NekFFTW", m_homogeneousBasis->GetNumPoints());
98
            }
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

            if(m_dealiasing)
            {
                if(m_useFFT)
                {
                    NekDouble size = 1.5*m_homogeneousBasis->GetNumPoints();
                    m_padsize = int(size);
                    m_FFT_deal = LibUtilities::GetNektarFFTFactory()
                                    .CreateInstance("NekFFTW", m_padsize);
                }
                else
                {
                    ASSERTL0(false, "Dealiasing available just in combination "
                                    "with FFTW");
                }
            }
        }
116
117
118
119
120
121
122


        /**
         * @param   In          ExpListHomogeneous1D object to copy.
         */
        ExpListHomogeneous1D::ExpListHomogeneous1D(const ExpListHomogeneous1D &In):
            ExpList(In,false),
123
            m_transposition(In.m_transposition),
124
            m_StripZcomm(In.m_StripZcomm),
125
            m_useFFT(In.m_useFFT),
126
            m_FFT(In.m_FFT),
127
            m_tmpIN(In.m_tmpIN),
128
            m_tmpOUT(In.m_tmpOUT),
129
130
131
132
            m_homogeneousBasis(In.m_homogeneousBasis),
            m_lhom(In.m_lhom), 
            m_homogeneous1DBlockMat(In.m_homogeneous1DBlockMat),
            m_dealiasing(In.m_dealiasing),
133
            m_padsize(In.m_padsize)
134
135
136
        {
            m_planes = Array<OneD, ExpListSharedPtr>(In.m_planes.num_elements());
        }
137
138
139
140
141
142
143
        
        ExpListHomogeneous1D::ExpListHomogeneous1D(const ExpListHomogeneous1D &In,
                                            const std::vector<unsigned int> &eIDs):
            ExpList(In,eIDs,false),
            m_transposition(In.m_transposition),
            m_useFFT(In.m_useFFT),
            m_FFT(In.m_FFT),
144
            m_tmpIN(In.m_tmpIN),
145
            m_tmpOUT(In.m_tmpOUT),
146
147
            m_homogeneousBasis(In.m_homogeneousBasis),
            m_lhom(In.m_lhom), 
148
            m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr()),
149
            m_dealiasing(In.m_dealiasing),
150
            m_padsize(In.m_padsize)
151
152
153
154
155
156
157
158
159
160
        {
            m_planes = Array<OneD, ExpListSharedPtr>(In.m_planes.num_elements());
        }

        /**
         * Destructor
         */
        ExpListHomogeneous1D::~ExpListHomogeneous1D()
        {
        }
Yan Bao's avatar
Yan Bao committed
161
    
162
        void ExpListHomogeneous1D::v_HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray, 
163
164
165
166
                                                         Array<OneD, NekDouble> &outarray, 
                                                         CoeffState coeffstate,
                                                         bool Shuff,
                                                         bool UnShuff)
167
        {
168
            // Forwards trans
169
            Homogeneous1DTrans(inarray,outarray,true,coeffstate,Shuff,UnShuff);
170
        }
Yan Bao's avatar
Yan Bao committed
171
    
172
        void ExpListHomogeneous1D::v_HomogeneousBwdTrans(const Array<OneD, const NekDouble> &inarray, 
173
174
175
176
                                                         Array<OneD, NekDouble> &outarray, 
                                                         CoeffState coeffstate,
                                                         bool Shuff,
                                                         bool UnShuff)
177
        {
178
            // Backwards trans
179
            Homogeneous1DTrans(inarray,outarray,false,coeffstate,Shuff,UnShuff);
180
        }
Yan Bao's avatar
Yan Bao committed
181
        
182
183
184
185
186
187
188
189
        /**
         * Dealiasing routine
         */
        void ExpListHomogeneous1D::v_DealiasedProd(const Array<OneD, NekDouble> &inarray1,
                                                   const Array<OneD, NekDouble> &inarray2,
                                                   Array<OneD, NekDouble> &outarray, 
                                                   CoeffState coeffstate)
        {
190
191
            // inarray1 = first term of the product in full physical space
            // inarray2 = second term of the product in full physical space
192
            // dealiased product stored in outarray
193
194
195
196
197

            int num_dofs = inarray1.num_elements();

            int N = m_homogeneousBasis->GetNumPoints();

198
199
            Array<OneD, NekDouble> V1(num_dofs);
            Array<OneD, NekDouble> V2(num_dofs);
200
201
202
203
204
205
            Array<OneD, NekDouble> V1V2(num_dofs);

            HomogeneousFwdTrans(inarray1,V1,coeffstate);
            HomogeneousFwdTrans(inarray2,V2,coeffstate);

            int num_points_per_plane = num_dofs/m_planes.num_elements();
Yan Bao's avatar
Yan Bao committed
206
207
208
209
210
211
212
213
214
            int num_proc;
            if(!m_session->DefinesSolverInfo("HomoStrip"))
            {
                num_proc             = m_comm->GetColumnComm()->GetSize();
            }
            else
            {
                num_proc             = m_StripZcomm->GetSize();
            }
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
            int num_dfts_per_proc    = num_points_per_plane / num_proc
                                        + (num_points_per_plane % num_proc > 0);

            Array<OneD, NekDouble> ShufV1(num_dfts_per_proc*N,0.0);
            Array<OneD, NekDouble> ShufV2(num_dfts_per_proc*N,0.0);
            Array<OneD, NekDouble> ShufV1V2(num_dfts_per_proc*N,0.0);

            Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize,0.0);
            Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize,0.0);
            Array<OneD, NekDouble> ShufV1_PAD_phys(m_padsize,0.0);
            Array<OneD, NekDouble> ShufV2_PAD_phys(m_padsize,0.0);

            Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize,0.0);
            Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize,0.0);

            m_transposition->Transpose(V1, ShufV1, false, LibUtilities::eXYtoZ);
            m_transposition->Transpose(V2, ShufV2, false, LibUtilities::eXYtoZ);

            // Looping on the pencils
234
            for(int i = 0 ; i < num_dfts_per_proc ; i++)
235
            {
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
                // Copying the i-th pencil pf lenght N into a bigger
                // pencil of lenght 2N We are in Fourier space
                Vmath::Vcopy(N, &(ShufV1[i*N]), 1, &(ShufV1_PAD_coef[0]), 1);
                Vmath::Vcopy(N, &(ShufV2[i*N]), 1, &(ShufV2_PAD_coef[0]), 1);

                // Moving to physical space using the padded system
                m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys);
                m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys);

                // Perfroming the vectors multiplication in physical space on
                // the padded system
                Vmath::Vmul(m_padsize, ShufV1_PAD_phys,   1,
                                       ShufV2_PAD_phys,   1,
                                       ShufV1V2_PAD_phys, 1);

                // Moving back the result (V1*V2)_phys in Fourier space, padded
                // system
                m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);

                // Copying the first half of the padded pencil in the full
                // vector (Fourier space)
                Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1,
                                &(ShufV1V2[i*N]),        1);
259
            }
260
261
262
263
264
265

            m_transposition->Transpose(ShufV1V2, V1V2, false,
                                       LibUtilities::eZtoXY);

            // Moving the results in physical space for the output
            HomogeneousBwdTrans(V1V2, outarray, coeffstate);
266
        }
267

268
269
270
271
        /**
         * Forward transform
         */
        void ExpListHomogeneous1D::v_FwdTrans(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray, CoeffState coeffstate )
272
273
274
        {
            int cnt = 0, cnt1 = 0;
            Array<OneD, NekDouble> tmparray;
275
            
276
            for(int n = 0; n < m_planes.num_elements(); ++n)
277
            {
278
                m_planes[n]->FwdTrans(inarray+cnt, tmparray = outarray + cnt1,
279
                                      coeffstate);
280
                cnt   += m_planes[n]->GetTotPoints();
281
282
                
                cnt1  += m_planes[n]->GetNcoeffs(); // need to skip ncoeffs
283
            }
284
            if(!m_WaveSpace)
285
286
287
            {
                HomogeneousFwdTrans(outarray,outarray,coeffstate);
            }
288
289
        }

290
291
292
        /**
         * Forward transform element by element
         */
293
294
295
296
        void ExpListHomogeneous1D::v_FwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
        {
            int cnt = 0, cnt1 = 0;
            Array<OneD, NekDouble> tmparray;
297
            
298
            //spectral element FwdTrans plane by plane
299
            for(int n = 0; n < m_planes.num_elements(); ++n)
300
301
302
303
304
            {
                m_planes[n]->FwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
                cnt   += m_planes[n]->GetTotPoints();
                cnt1  += m_planes[n]->GetNcoeffs();
            }
305
306
307
308
            if(!m_WaveSpace)
            {
                HomogeneousFwdTrans(outarray,outarray);
            }
309
        }
310
311
312
313
314
        
        /**
         * Backward transform
         */
        void ExpListHomogeneous1D::v_BwdTrans(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray, CoeffState coeffstate)
315
316
317
        {
            int cnt = 0, cnt1 = 0;
            Array<OneD, NekDouble> tmparray;
318
            
319
            for(int n = 0; n < m_planes.num_elements(); ++n)
320
321
            {
                m_planes[n]->BwdTrans(inarray+cnt, tmparray = outarray + cnt1,
322
323
324
325
326
327
328
                                      coeffstate);
                cnt  += m_planes[n]->GetNcoeffs();
                cnt1 += m_planes[n]->GetTotPoints();
            }
            if(!m_WaveSpace)
            {
                HomogeneousBwdTrans(outarray,outarray);
329
330
            }
        }
Yan Bao's avatar
Yan Bao committed
331
    
332
333
334
335
        /**
         * Backward transform element by element
         */
        void ExpListHomogeneous1D::v_BwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
336
337
338
        {
            int cnt = 0, cnt1 = 0;
            Array<OneD, NekDouble> tmparray;
339
            
340
            for(int n = 0; n < m_planes.num_elements(); ++n)
341
342
            {
                m_planes[n]->BwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
343
344
                
                cnt    += m_planes[n]->GetNcoeffs();
345
346
                cnt1   += m_planes[n]->GetTotPoints();
            }
347
348
349
350
            if(!m_WaveSpace)
            {
                HomogeneousBwdTrans(outarray,outarray);
            }
351
        }
352
353
354
355
356
        
        /**
         * Inner product
         */
        void ExpListHomogeneous1D::v_IProductWRTBase(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray, CoeffState coeffstate)
357
358
359
        {
            int cnt = 0, cnt1 = 0;
            Array<OneD, NekDouble> tmparray;
360
            
361
            for(int n = 0; n < m_planes.num_elements(); ++n)
362
            {
363
                m_planes[n]->IProductWRTBase(inarray+cnt, tmparray = outarray + cnt1,coeffstate);
364

365
                cnt1    += m_planes[n]->GetNcoeffs();
366
                cnt   += m_planes[n]->GetTotPoints();
367
368
            }
        }
Yan Bao's avatar
Yan Bao committed
369
    
370
371
372
373
        /**
         * Inner product element by element
         */
        void ExpListHomogeneous1D::v_IProductWRTBase_IterPerExp(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
374
        { 
375
376
            int cnt = 0, cnt1 = 0;
            Array<OneD, NekDouble> tmparray;
377
            
378
            for(int n = 0; n < m_planes.num_elements(); ++n)
379
380
            {
                m_planes[n]->IProductWRTBase_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
Yan Bao's avatar
Yan Bao committed
381
        
382
383
                cnt1  += m_planes[n]->GetNcoeffs();
                cnt   += m_planes[n]->GetTotPoints();
384
            } 
385
        }
Yan Bao's avatar
Yan Bao committed
386
    
387
388
389
390
391
392
393
394
        /**
         * Homogeneous transform Bwd/Fwd (MVM and FFT)
         */
        void ExpListHomogeneous1D::Homogeneous1DTrans(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray, 
                                                      bool IsForwards, 
                                                      CoeffState coeffstate,
                                                      bool Shuff,
                                                      bool UnShuff)
395
        {
396
397
            int num_dofs;
            
398
            if(IsForwards) 
399
400
401
402
403
404
405
406
            {
                num_dofs = inarray.num_elements();
            }
            else
            {
                num_dofs = outarray.num_elements();
            }
            
Alessandro Bolis's avatar
fftw    
Alessandro Bolis committed
407
            if(m_useFFT)
Yan Bao's avatar
Yan Bao committed
408
            {        
409
                int num_points_per_plane = num_dofs/m_planes.num_elements();
Yan Bao's avatar
Yan Bao committed
410
                int num_dfts_per_proc;
Yan Bao's avatar
Yan Bao committed
411
412
413
414
415
416
417
418
419
420
421
422
                if(!m_session->DefinesSolverInfo("HomoStrip"))
                {
                    int nP = m_comm->GetColumnComm()->GetSize();
                    num_dfts_per_proc = num_points_per_plane / nP
                                      + (num_points_per_plane % nP > 0);
                }
                else
                {
                    int nP = m_StripZcomm->GetSize();
                    num_dfts_per_proc = num_points_per_plane / nP
                                      + (num_points_per_plane % nP > 0);
                }
423

424
                Array<OneD, NekDouble> fft_in (num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),0.0);
425
                Array<OneD, NekDouble> fft_out(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),0.0);
Yan Bao's avatar
Yan Bao committed
426
        
427
428
429
430
431
432
                if(Shuff)
                {
                    m_transposition->Transpose(inarray,fft_in,false,LibUtilities::eXYtoZ);
                }
                else 
                {
433
434
                    Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
                                 inarray,1,fft_in,1);
435
436
                }
                
437
438
                if(IsForwards)
                {
439
                    for(int i = 0 ; i < num_dfts_per_proc ; i++)
440
                    {
441
                        m_FFT->FFTFwdTrans(m_tmpIN = fft_in + i*m_homogeneousBasis->GetNumPoints(), m_tmpOUT = fft_out + i*m_homogeneousBasis->GetNumPoints());
442
443
444
445
                    }
                }
                else 
                {
446
                    for(int i = 0 ; i < num_dfts_per_proc ; i++)
447
                    {
448
                        m_FFT->FFTBwdTrans(m_tmpIN = fft_in + i*m_homogeneousBasis->GetNumPoints(), m_tmpOUT = fft_out + i*m_homogeneousBasis->GetNumPoints());
449
450
                    }
                }
Yan Bao's avatar
Yan Bao committed
451
        
452
453
454
455
456
457
                if(UnShuff)
                {
                    m_transposition->Transpose(fft_out,outarray,false,LibUtilities::eZtoXY);
                }
                else 
                {
458
459
                    Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
                                 fft_out,1,outarray,1);
460
                }
461
462
463
464
            }
            else 
            {
                DNekBlkMatSharedPtr blkmat;
Yan Bao's avatar
Yan Bao committed
465
        
466
                if(num_dofs == m_npoints) //transform phys space
467
                {
468
469
470
471
472
473
474
475
                    if(IsForwards)
                    {
                        blkmat = GetHomogeneous1DBlockMatrix(eForwardsPhysSpace1D);
                    }
                    else
                    {
                        blkmat = GetHomogeneous1DBlockMatrix(eBackwardsPhysSpace1D);
                    }
476
                }
477
                else
478
                {
479
480
                    if(IsForwards)
                    {
481
                        blkmat = GetHomogeneous1DBlockMatrix(eForwardsCoeffSpace1D,coeffstate);
482
483
484
                    }
                    else
                    {
485
                        blkmat = GetHomogeneous1DBlockMatrix(eBackwardsCoeffSpace1D,coeffstate);
486
                    }
487
                }
Yan Bao's avatar
Yan Bao committed
488
        
Alessandro Bolis's avatar
fftw    
Alessandro Bolis committed
489
490
                int nrows = blkmat->GetRows();
                int ncols = blkmat->GetColumns();
Yan Bao's avatar
Yan Bao committed
491
        
492
493
                Array<OneD, NekDouble> sortedinarray(ncols,0.0);
                Array<OneD, NekDouble> sortedoutarray(nrows,0.0);
Yan Bao's avatar
Yan Bao committed
494
        
495
496
497
498
499
500
501
502
503
                if(Shuff)
                {
                    m_transposition->Transpose(inarray,sortedinarray,!IsForwards,LibUtilities::eXYtoZ);
                }
                else 
                {
                    Vmath::Vcopy(ncols,inarray,1,sortedinarray,1);
                }
                
Alessandro Bolis's avatar
fftw    
Alessandro Bolis committed
504
                // Create NekVectors from the given data arrays
505
                NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
506
                NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
Yan Bao's avatar
Yan Bao committed
507
        
Alessandro Bolis's avatar
fftw    
Alessandro Bolis committed
508
509
                // Perform matrix-vector multiply.
                out = (*blkmat)*in;
Yan Bao's avatar
Yan Bao committed
510
        
511
512
513
514
515
516
                if(UnShuff)
                {
                    m_transposition->Transpose(sortedoutarray,outarray,IsForwards,LibUtilities::eZtoXY);
                }
                else 
                {
Douglas Serson's avatar
Douglas Serson committed
517
                    Vmath::Vcopy(nrows,sortedoutarray,1,outarray,1);
518
519
                }
                
520
            }
521
522
        }

523
        DNekBlkMatSharedPtr ExpListHomogeneous1D::GetHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate) const
524
525
        {
            Homo1DBlockMatrixMap::iterator matrixIter = m_homogeneous1DBlockMat->find(mattype);
526
            
527
528
            if(matrixIter == m_homogeneous1DBlockMat->end())
            {
529
                return ((*m_homogeneous1DBlockMat)[mattype] =
530
                        GenHomogeneous1DBlockMatrix(mattype,coeffstate));
531
532
533
534
535
536
537
            }
            else
            {
                return matrixIter->second;
            }
        }

538

539
        DNekBlkMatSharedPtr ExpListHomogeneous1D::GenHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate) const
540
        {
541
542
            DNekMatSharedPtr    loc_mat;
            DNekBlkMatSharedPtr BlkMatrix;
543
544
545
            int n_exp = 0;
            int num_trans_per_proc = 0;
            
546
547
            if((mattype == eForwardsCoeffSpace1D)
               ||(mattype == eBackwardsCoeffSpace1D)) // will operate on m_coeffs
548
            {
549
                n_exp = m_planes[0]->GetNcoeffs();
550
551
552
553
554
            }
            else
            {
                n_exp = m_planes[0]->GetTotPoints(); // will operatore on m_phys
            }
555

556
            num_trans_per_proc = n_exp/m_comm->GetColumnComm()->GetSize() + (n_exp%m_comm->GetColumnComm()->GetSize() > 0);
557

558
559
            Array<OneD,unsigned int> nrows(num_trans_per_proc);
            Array<OneD,unsigned int> ncols(num_trans_per_proc);
560

561
            if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
562
            {
563
564
                nrows = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumModes());
                ncols = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumPoints());
565
566
567
            }
            else
            {
568
569
                nrows = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumPoints());
                ncols = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumModes());
570
571
572
573
574
            }

            MatrixStorage blkmatStorage = eDIAGONAL;
            BlkMatrix = MemoryManager<DNekBlkMat>
                ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
575

Yan Bao's avatar
Yan Bao committed
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
            //Half Mode
            if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
            {
                StdRegions::StdPointExp StdPoint(m_homogeneousBasis->GetBasisKey());
                
                if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
                {
                    StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
                                                    StdPoint.DetShapeType(),
                                                    StdPoint);
                    
                    loc_mat = StdPoint.GetStdMatrix(matkey);
                }
                else
                {
                    StdRegions::StdMatrixKey matkey(StdRegions::eBwdTrans,
                                                    StdPoint.DetShapeType(),
                                                    StdPoint);
                    
                    loc_mat = StdPoint.GetStdMatrix(matkey);
                }
            }
            //other cases
            else 
            {
                StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
                
                if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
                {
                    StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
                                                    StdSeg.DetShapeType(),
                                                    StdSeg);
                    
                    loc_mat = StdSeg.GetStdMatrix(matkey);
                }
                else
                {
                    StdRegions::StdMatrixKey matkey(StdRegions::eBwdTrans,
                                                    StdSeg.DetShapeType(),
                                                    StdSeg);
                    
                    loc_mat = StdSeg.GetStdMatrix(matkey);
                }                
                
            }
621

622
            // set up array of block matrices.
623
            for(int i = 0; i < num_trans_per_proc; ++i)
624
625
626
            {
                BlkMatrix->SetBlock(i,i,loc_mat);
            }
627

628
629
630
            return BlkMatrix;
        }

631
        std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpListHomogeneous1D::v_GetFieldDefinitions()
632
        {
633
            std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
634
            
Yan Bao's avatar
Yan Bao committed
635
            // Set up Homogeneous length details.
636
            Array<OneD,LibUtilities::BasisSharedPtr> HomoBasis(1,m_homogeneousBasis);
637
            
638
            std::vector<NekDouble> HomoLen;
639
            HomoLen.push_back(m_lhom);
640
            
641
642
            std::vector<unsigned int> StripsIDs;

643
644
645
646
647
648
            bool strips;
            m_session->MatchSolverInfo("HomoStrip","True",strips,false);
            if (strips)
            {
                StripsIDs.push_back(m_transposition->GetStripID());
            }
649
            
650
            std::vector<unsigned int> PlanesIDs;
651
652
653
            int IDoffset = 0;

            // introduce a 2 plane offset for single mode case so can
654
            // be post-processed or used in MultiMode expansion.
655
656
657
658
659
            if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode)
            {
                IDoffset  = 2;
            }

660
661
            for(int i = 0; i < m_planes.num_elements(); i++)
            {
662
                PlanesIDs.push_back(m_transposition->GetPlaneID(i)+IDoffset);
663
664
            }
            
665
666
            m_planes[0]->GeneralGetFieldDefinitions(returnval, 1, HomoBasis, 
                    HomoLen, strips, StripsIDs, PlanesIDs);
667
            return returnval;
668
669
        }

670
        void  ExpListHomogeneous1D::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
671
        {
672
            // Set up Homogeneous length details.
673
            Array<OneD,LibUtilities::BasisSharedPtr> HomoBasis(1,m_homogeneousBasis);
674
            
675
            std::vector<NekDouble> HomoLen;
676
            HomoLen.push_back(m_lhom);
677
            
678
679
            std::vector<unsigned int> StripsIDs;

680
681
682
683
684
685
            bool strips;
            m_session->MatchSolverInfo("HomoStrip","True",strips,false);
            if (strips)
            {
                StripsIDs.push_back(m_transposition->GetStripID());
            }
686
            
687
            std::vector<unsigned int> PlanesIDs;
688
689
690
691
692
693
694
            int IDoffset = 0;

            if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode)
            {
                IDoffset = 2;
            }

695
696
            for(int i = 0; i < m_planes.num_elements(); i++)
            {
697
                PlanesIDs.push_back(m_transposition->GetPlaneID(i)+IDoffset);
698
699
700
            }
            
            // enforce NumHomoDir == 1 by direct call
701
702
            m_planes[0]->GeneralGetFieldDefinitions(fielddef, 1, HomoBasis,
                    HomoLen, strips, StripsIDs, PlanesIDs);
703
704
705
        }


706
707
708
709
710
711
        /** This routine appends the data from the expansion list into
            the output format where each element is given by looping
            over its Fourier modes where as data in the expandion is
            stored with all consecutive elements and then the Fourier
            modes
         */ 
712
        void ExpListHomogeneous1D::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
713
714
715
        {
            int i,n;
            int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
716

717
            // Determine mapping from element ids to location in
718
            // expansion list
719
            if (m_elmtToExpId.size() == 0)
720
            {
721
722
723
724
                for(i = 0; i < m_planes[0]->GetExpSize(); ++i)
                {
                    m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
                }
725
            }
726

727
            for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
728
            {
729
                int eid     = m_elmtToExpId[fielddef->m_elementIDs[i]];
730
731
                int datalen = (*m_exp)[eid]->GetNcoeffs();

732
                for(n = 0; n < m_planes.num_elements(); ++n)
733
                {
734
                    fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]+n*ncoeffs_per_plane],&coeffs[m_coeff_offset[eid]+n*ncoeffs_per_plane]+datalen);
735
736
737
                }
            }
        }
Yan Bao's avatar
Yan Bao committed
738
        
739
        void ExpListHomogeneous1D::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
740
741
742
        {
           v_AppendFieldData(fielddef,fielddata,m_coeffs);
        }
743
744

        //Extract the data in fielddata into the m_coeff list
745
        void ExpListHomogeneous1D::v_ExtractDataToCoeffs(
746
747
748
749
            LibUtilities::FieldDefinitionsSharedPtr    &fielddef,
            std::vector<NekDouble>       &fielddata,
            std::string                  &field,
            Array<OneD, NekDouble>       &coeffs)
750
751
        {
            int i,n;
752
753
            int offset  = 0;
            int nzmodes = 1;
754
            int datalen = fielddata.size()/fielddef->m_fields.size();
755
            std::vector<unsigned int> fieldDefHomoZids;
756
757
            
            
758
            // Find data location according to field definition
759
            for(i = 0; i < fielddef->m_fields.size(); ++i)
760
            {
761
                if(fielddef->m_fields[i] == field)
762
763
764
765
766
                {
                    break;
                }
                offset += datalen;
            }
767

768
769
770
            if(i == fielddef->m_fields.size())
            {
                cout << "Field "<< field<< "not found in data file. "  << endl;
771
772
            }
            else
773
            {
774
                
775
776
777
                int modes_offset = 0;
                int planes_offset = 0;
                Array<OneD, NekDouble> coeff_tmp;
778
                boost::unordered_map<int,int>::iterator it;
779
                
780
781
782
                // Build map of plane IDs lying on this processor and determine
                // mapping from element ids to location in expansion list.
                if (m_zIdToPlane.size() == 0)
783
                {
784
785
786
787
788
789
790
791
792
                    for (i = 0; i < m_planes.num_elements(); ++i)
                    {
                        m_zIdToPlane[m_transposition->GetPlaneID(i)] = i;
                    }

                    for (i = 0; i < m_planes[0]->GetExpSize(); ++i)
                    {
                        m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
                    }
793
                }
794
                
795
                if(fielddef->m_numHomogeneousDir)
796
                {
797
                    nzmodes = fielddef->m_homogeneousZIDs.size();
798
799
800
801
802
803
804
805
806
807
808
809
810
                    fieldDefHomoZids = fielddef->m_homogeneousZIDs;
                }
                else // input file is 2D and so set nzmodes to 1
                {
                    nzmodes = 1;
                    fieldDefHomoZids.push_back(0);
                }

                // calculate number of modes in the current partition
                int ncoeffs_per_plane = m_planes[0]->GetNcoeffs(); 
                
                for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
                {
811
                    if(fielddef->m_uniOrder == true) // reset modes_offset to zero
812
                    {
813
                        modes_offset = 0;
814
815
                    }
                    
816
817
818
819
                    int datalen = LibUtilities::GetNumberOfCoefficients(fielddef->m_shapeType, 
                                                                        fielddef->m_numModes,
                                                                        modes_offset);

820
                    it = m_elmtToExpId.find(fielddef->m_elementIDs[i]);
821
822
                    
                    // ensure element is on this partition for parallel case. 
823
                    if(it == m_elmtToExpId.end())
824
825
826
                    {
                        // increase offset for correct FieldData access
                        offset += datalen*nzmodes;
827
828
                        modes_offset += (*m_exp)[0]->GetNumBases() +
                                        fielddef->m_numHomogeneousDir;
829
830
831
832
833
834
                        continue;
                    }
                    
                    int eid = it->second;
                        
                    
835
                    for(n = 0; n < nzmodes; ++n, offset += datalen)
836
                    {
837
                        
838
                        it = m_zIdToPlane.find(fieldDefHomoZids[n]);
839
                        
840
                        // Check to make sure this mode number lies in this field.
841
                        if (it == m_zIdToPlane.end())
842
843
                        {
                            continue;
844
845
                        } 
                        
846
847
848
849
                        planes_offset = it->second;
                        if(datalen == (*m_exp)[eid]->GetNcoeffs())
                        {
                            Vmath::Vcopy(datalen,&fielddata[offset],1,&coeffs[m_coeff_offset[eid]+planes_offset*ncoeffs_per_plane],1);
850
                        }
851
852
                        else // unpack data to new order
                        {
853
                            (*m_exp)[eid]->ExtractDataToCoeffs(&fielddata[offset], fielddef->m_numModes,modes_offset,&coeffs[m_coeff_offset[eid] + planes_offset*ncoeffs_per_plane]);
854
                        }
855
                    }
856
                    modes_offset += (*m_exp)[0]->GetNumBases() + fielddef->m_numHomogeneousDir;
857
858
859
                }
            }
        }
Yan Bao's avatar
Yan Bao committed
860
        
861
862
863
        //Extract the data in fielddata into the m_coeff list
        void ExpListHomogeneous1D::v_ExtractCoeffsToCoeffs(
                                                           const boost::shared_ptr<ExpList> &fromExpList,const  Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
864
        {
865
            int i;
866
            int fromNcoeffs_per_plane = fromExpList->GetPlane(0)->GetNcoeffs();
867
            int toNcoeffs_per_plane = m_planes[0]->GetNcoeffs();
868
869
870
            Array<OneD, NekDouble> tocoeffs_tmp, fromcoeffs_tmp; 
            
            for(i = 0; i < m_planes.num_elements(); ++i)
871
            {
872
                m_planes[i]->ExtractCoeffsToCoeffs(fromExpList->GetPlane(i),fromcoeffs_tmp =  fromCoeffs + fromNcoeffs_per_plane*i, tocoeffs_tmp = toCoeffs + toNcoeffs_per_plane*i);
873
874
            }
        }
875

876
        void ExpListHomogeneous1D::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
877
878
                                        std::string var)
        {
879
880
881
882
883
884
885
            // If there is only one plane (e.g. HalfMode), we write a 2D plane.
            if (m_planes.num_elements() == 1)
            {
                m_planes[0]->WriteVtkPieceData(outfile, expansion, var);
                return;
            }

886
887
888
889
            int i;
            int nq = (*m_exp)[expansion]->GetTotPoints();
            int npoints_per_plane = m_planes[0]->GetTotPoints();

890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
            // If we are using Fourier points, output extra plane to fill domain
            int outputExtraPlane = 0;
            Array<OneD, NekDouble> extraPlane;
            if ( m_homogeneousBasis->GetBasisType()   == LibUtilities::eFourier
               && m_homogeneousBasis->GetPointsType() ==
                    LibUtilities::eFourierEvenlySpaced)
            {
                outputExtraPlane = 1;
                // Get extra plane data
                if (m_StripZcomm->GetSize() == 1)
                {
                    extraPlane = m_phys + m_phys_offset[expansion];
                }
                else
                {
                    // Determine to and from rank for communication
                    int size     = m_StripZcomm->GetSize();
                    int rank     = m_StripZcomm->GetRank();
                    int fromRank = (rank+1) % size;
                    int toRank   = (rank == 0) ? size-1 : rank-1;
                    // Communicate using SendRecv
                    extraPlane = Array<OneD, NekDouble>(nq);
                    Array<OneD, NekDouble> send (nq,
                            m_phys + m_phys_offset[expansion]);
                    m_StripZcomm->SendRecv(toRank, send,
                                           fromRank, extraPlane);
                }
            }

919
            // printing the fields of that zone
920
            outfile << "        <DataArray type=\"Float64\" Name=\""
921
922
923
924
925
926
927
928
929
930
                    << var << "\">" << endl;
            outfile << "          ";
            for (int n = 0; n < m_planes.num_elements(); ++n)
            {
                const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion] + n*npoints_per_plane;
                for(i = 0; i < nq; ++i)
                {
                    outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
                }
            }
931
932
933
934
935
936
937
938
            if (outputExtraPlane)
            {
                for(i = 0; i < nq; ++i)
                {
                    outfile << (fabs(extraPlane[i]) < NekConstants::kNekZeroTol ?
                                0 : extraPlane[i]) << " ";
                }
            }
939
940
941
            outfile << endl;
            outfile << "        </DataArray>" << endl;
        }
942
        
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
        void ExpListHomogeneous1D::v_PhysInterp1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
        {
            int cnt,cnt1;
            Array<OneD, NekDouble> tmparray;
            cnt  = m_planes[0]->GetTotPoints();
            cnt1 = m_planes[0]->Get1DScaledTotPoints(scale);
            
            ASSERTL1(m_planes.num_elements()*cnt1 <= outarray.num_elements(),"size of outarray does not match internal estimage");
            
            
            for(int i = 0; i < m_planes.num_elements(); i++)
            {
         
                m_planes[i]->PhysInterp1DScaled(scale,inarray+i*cnt,
                                                 tmparray = outarray+i*cnt1);
            }
        }


        void ExpListHomogeneous1D::v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
        {
            int cnt,cnt1;
            Array<OneD, NekDouble> tmparray;
            cnt  = m_planes[0]->Get1DScaledTotPoints(scale);
            cnt1 = m_planes[0]->GetTotPoints();
            
            ASSERTL1(m_planes.num_elements()*cnt <= inarray.num_elements(),"size of outarray does not match internal estimage");
            
            
            for(int i = 0; i < m_planes.num_elements(); i++)
            {
                m_planes[i]->PhysGalerkinProjection1DScaled(scale,inarray+i*cnt,
                                                 tmparray = outarray+i*cnt1);
            }
            
        }
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
        void ExpListHomogeneous1D::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
                                               Array<OneD, NekDouble> &out_d0,
                                               Array<OneD, NekDouble> &out_d1, 
                                               Array<OneD, NekDouble> &out_d2)
        {
            int nT_pts = inarray.num_elements();          //number of total points = n. of Fourier points * n. of points per plane (nT_pts)
            int nP_pts = nT_pts/m_planes.num_elements();    //number of points per plane = n of Fourier transform required (nP_pts)
            
            Array<OneD, NekDouble> temparray(nT_pts);
            Array<OneD, NekDouble> outarray(nT_pts);
            Array<OneD, NekDouble> tmp1;
            Array<OneD, NekDouble> tmp2;
            Array<OneD, NekDouble> tmp3;            
            
            for(int i = 0; i < m_planes.num_elements(); i++)
            {
                m_planes[i]->PhysDeriv(inarray + i*nP_pts ,tmp2 = out_d0 + i*nP_pts , tmp3 = out_d1 + i*nP_pts );
            }
            
998
            if(out_d2 != NullNekDouble1DArray)
999
            {
1000
                if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode || 
Yan Bao's avatar
Yan Bao committed
1001
                   m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)            
1002
                {
1003
1004
1005
1006
1007
1008
1009
1010
                    if(m_WaveSpace)
                    {
                        temparray = inarray;
                    }
                    else 
                    { 
                        HomogeneousFwdTrans(inarray,temparray);
                    }
1011
                    
1012
1013
                    NekDouble sign = -1.0;
                    NekDouble beta;
1014
                    
1015
1016
                    //Half Mode
                    if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe)
1017
                    {
1018
1019
1020
1021
1022
1023
1024
1025
1026
                        beta = sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
                        
                        Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
                    }
                    else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
                    {
                        beta = -sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
                        
                        Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
1027
                    }
1028
                    
1029
1030
                    //Fully complex
                    else
1031
                    {
1032
1033
1034
1035
1036
1037
1038
1039
                        for(int i = 0; i < m_planes.num_elements(); i++)
                        {
                            beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
                            
                            Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
                            
                            sign = -1.0*sign;
                        }