DriverSteadyState.cpp 12.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
///////////////////////////////////////////////////////////////////////////////
//
// File DriverSteadyState.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: Incompressible Navier Stokes solver
//
///////////////////////////////////////////////////////////////////////////////

#include <SolverUtils/DriverSteadyState.h>
Bastien Jordi's avatar
Bastien Jordi committed
37
#include <SolverUtils/AdvectionSystem.h>
38
39
40
41


namespace Nektar
{
Bastien Jordi's avatar
Bastien Jordi committed
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
    namespace SolverUtils
    {
        string DriverSteadyState::className = GetDriverFactory().RegisterCreatorFunction("SteadyState", DriverSteadyState::create);
        string DriverSteadyState::driverLookupId = LibUtilities::SessionReader::RegisterEnumValue("Driver","SteadyState",0);
        
        /**
         * 
         */
        DriverSteadyState::DriverSteadyState(const LibUtilities::SessionReaderSharedPtr pSession)
        : Driver(pSession)
        {
        }
        
        
        /**
         * 
         */
        DriverSteadyState:: ~DriverSteadyState()
        {
        }
        
        
        /**
         * 
         */
        void DriverSteadyState::v_InitObject(ostream &out)
        {
            Driver::v_InitObject(out);
        }
        
        
        void DriverSteadyState::v_Execute(ostream &out)
        
        {
            //With a loop over "DoSolve", this Driver implements the "encaplulated" Selective Frequency Damping method
            //to find the steady state of a flow above the critical Reynolds number.
            
            m_equ[0]->PrintSummary(out);
            m_equ[0]->DoInitialise();
            
Bastien Jordi's avatar
Bastien Jordi committed
82
83
84
85
86
87
88
            
            ////////////////////////////////////////////////////////////////////////////////////////////////////////////
            AdvectionSystemSharedPtr A = boost::dynamic_pointer_cast<AdvectionSystem>(m_equ[0]);
            ////////////////////////////////////////////////////////////////////////////////////////////////////////////
            
            
            
89
            // - SFD Routine -
90
91
92
93
            // Compressible case
            NumVar_SFD = m_equ[0]->UpdateFields()[0]->GetCoordim(0);            
            if (m_session->GetSolverInfo("EqType") == "EulerCFE" || 
                m_session->GetSolverInfo("EqType") == "NavierStokesCFE")
Daniele de Grazia's avatar
Daniele de Grazia committed
94
            {
95
                NumVar_SFD += 2; //Number of variables for the compressible equations
Daniele de Grazia's avatar
Daniele de Grazia committed
96
97
98
99
100
101
            }
                
            Array<OneD, Array<OneD, NekDouble> > q0(NumVar_SFD);
            Array<OneD, Array<OneD, NekDouble> > q1(NumVar_SFD);
            Array<OneD, Array<OneD, NekDouble> > qBar0(NumVar_SFD);
            Array<OneD, Array<OneD, NekDouble> > qBar1(NumVar_SFD);
Bastien Jordi's avatar
Bastien Jordi committed
102
103
104
105
106
            
            NekDouble TOL(0);
            m_n=0;
            m_Check=0;			
            
107
            m_session->LoadParameter("FilterWidth", m_Delta0, 1);
Bastien Jordi's avatar
Bastien Jordi committed
108
            m_session->LoadParameter("ControlCoeff",m_X0, 1);
109
110
111
            m_session->LoadParameter("TOL", TOL, 1.0e-08);
            m_session->LoadParameter("IO_InfoSteps", m_infosteps, 1000);
            m_session->LoadParameter("IO_CheckSteps", m_checksteps, 100000);
Bastien Jordi's avatar
Bastien Jordi committed
112
113
            
            m_dt = m_equ[0]->GetTimeStep();
Bastien Jordi's avatar
Bastien Jordi committed
114
            
Bastien Jordi's avatar
Bastien Jordi committed
115
            //The time-step is applied here to clarify the matrix F expression
Bastien Jordi's avatar
Bastien Jordi committed
116
117
            m_X = m_X0*m_dt;
            m_Delta = m_Delta0/m_dt;
Bastien Jordi's avatar
Bastien Jordi committed
118
            
119
120
            //Exact solution of the Filters equation (ici on a X_tilde et Delta_tile!!!)
            c1 = 1.0/(1.0 + m_X*m_Delta);
Bastien Jordi's avatar
Bastien Jordi committed
121
122
123
124
125
            F11 = c1*(1.0 + m_X*m_Delta*exp(-(m_X + 1.0/m_Delta)));
            F12 = c1*(m_X*m_Delta*(1.0 - exp(-(m_X + 1.0/m_Delta))));
            F21 = c1*(1.0 - exp(-(m_X + 1.0/m_Delta)));
            F22 = c1*(m_X*m_Delta + exp(-(m_X + 1.0/m_Delta)));
            
Bastien Jordi's avatar
Bastien Jordi committed
126
            cout << "------------------ SFD Parameters ------------------" << endl;
Bastien Jordi's avatar
Bastien Jordi committed
127
128
            cout << "\tX = " << m_X0 << endl;
            cout << "\tDelta = " << m_Delta0 << endl;
Bastien Jordi's avatar
Bastien Jordi committed
129
130
            cout << "----------------------------------------------------" << endl;
            
131
            m_equ[0]->SetStepsToOne(); //m_steps is set to 1. Then "m_equ[0]->DoSolve()" will run for only one time step			
Bastien Jordi's avatar
Bastien Jordi committed
132
            ofstream m_file("ConvergenceHistory.txt", ios::out | ios::trunc);
Bastien Jordi's avatar
Bastien Jordi committed
133
            
Daniele de Grazia's avatar
Daniele de Grazia committed
134
            for(int i = 0; i < NumVar_SFD; ++i)
Bastien Jordi's avatar
Bastien Jordi committed
135
            {
Bastien Jordi's avatar
Bastien Jordi committed
136
                q0[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(), 0.0); //q0 is initialised
Bastien Jordi's avatar
Bastien Jordi committed
137
                qBar0[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(), 0.0);
Bastien Jordi's avatar
Bastien Jordi committed
138
                m_equ[0]->CopyFromPhysField(i, qBar0[i]); //qBar0 is initially set at beiing equal to the initial conditions provided in the input file
Bastien Jordi's avatar
Bastien Jordi committed
139
140
            }
            
141
            MaxNormDiff_q_qBar = 1.0;
142
            MaxNormDiff_q1_q0 = 1.0;
143
144
            Min_MaxNormDiff_q_qBar = MaxNormDiff_q_qBar;
            
Bastien Jordi's avatar
Bastien Jordi committed
145
            
Bastien Jordi's avatar
Bastien Jordi committed
146
            ////////////////////////////////////////////////////////////////////////////////////////////////////////////
147
            A->GetAdvObject()->SetBaseFlow(q0);
Bastien Jordi's avatar
Bastien Jordi committed
148
149
150
            ////////////////////////////////////////////////////////////////////////////////////////////////////////////
            
            
151
            while (max(MaxNormDiff_q_qBar, MaxNormDiff_q1_q0) > TOL)
Bastien Jordi's avatar
Bastien Jordi committed
152
            {
153
154
                //First order Splitting with exact resolution of the filters equation
                m_equ[0]->DoSolve();
155
                
Daniele de Grazia's avatar
Daniele de Grazia committed
156
                for(int i = 0; i < NumVar_SFD; ++i)
157
158
159
                {
                    m_equ[0]->CopyFromPhysField(i, q0[i]);
                    
Bastien Jordi's avatar
Bastien Jordi committed
160
                    //Apply the linear operator F to the outcome of the solver
161
162
163
164
165
                    ExactFilters(i, q0, qBar0, q1, qBar1);
                    
                    qBar0[i] = qBar1[i];
                    m_equ[0]->CopyToPhysField(i, q1[i]);     
                }
Bastien Jordi's avatar
Bastien Jordi committed
166
                
Bastien Jordi's avatar
Bastien Jordi committed
167
168
169
                
                if(m_infosteps && !((m_n+1)%m_infosteps))
                {
Bastien Jordi's avatar
Bastien Jordi committed
170
                    ConvergenceHistory(qBar1, q0, MaxNormDiff_q_qBar, MaxNormDiff_q1_q0);
Bastien Jordi's avatar
Bastien Jordi committed
171
172
173
174
175
176
177
178
                }
                
                if(m_checksteps && m_n&&(!((m_n+1)%m_checksteps)))
                {                      
                    m_Check++;                    
                    m_equ[0]->Checkpoint_Output(m_Check);                       
                }
                
Bastien Jordi's avatar
Bastien Jordi committed
179
180
181
                m_n++;
            }
            
Bastien Jordi's avatar
Bastien Jordi committed
182
183
184
            cout << "\nFINAL Filter Width: Delta = " << m_Delta << "; FINAL Control Coeff: X = " << m_X << "\n" << endl;
            m_Check++;
            m_equ[0]->Checkpoint_Output(m_Check);
Bastien Jordi's avatar
Bastien Jordi committed
185
            
Bastien Jordi's avatar
Bastien Jordi committed
186
187
188
189
190
191
192
193
194
195
            m_file.close();
            // - End SFD Routine -
            
            m_equ[0]->Output();
            
            // Evaluate and output computation time and solution accuracy.
            // The specific format of the error output is essential for the
            // regression tests to work.
            // Evaluate L2 Error
            for(int i = 0; i < m_equ[0]->GetNvariables(); ++i)
196
            {
Bastien Jordi's avatar
Bastien Jordi committed
197
198
199
200
201
202
203
                NekDouble vL2Error = m_equ[0]->L2Error(i,false);
                NekDouble vLinfError = m_equ[0]->LinfError(i);
                if (m_comm->GetRank() == 0)
                {
                    out << "L 2 error (variable " << m_equ[0]->GetVariable(i) << ") : " << vL2Error << endl;
                    out << "L inf error (variable " << m_equ[0]->GetVariable(i) << ") : " << vLinfError << endl;
                }
204
            }
Bastien Jordi's avatar
Bastien Jordi committed
205
206
207
208
209
210
211
212
213
214
215
        }
        
        
        void DriverSteadyState::ExactFilters(const int i,
                                             const Array<OneD, const Array<OneD, NekDouble> > &q0,
                                             const Array<OneD, const Array<OneD, NekDouble> > &qBar0,
                                             Array<OneD, Array<OneD, NekDouble> > &q1,
                                             Array<OneD, Array<OneD, NekDouble> > &qBar1)
        {
            q1[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(),0.0);
            qBar1[i] = Array<OneD, NekDouble> (m_equ[0]->GetTotPoints(),0.0);
216
            
Bastien Jordi's avatar
Bastien Jordi committed
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
            //Exact solution of the Filters equation            
            Vmath::Svtvp(q1[i].num_elements(), F11, q0[i], 1, q1[i], 1, q1[i], 1 );    
            Vmath::Svtvp(q1[i].num_elements(), F12, qBar0[i], 1, q1[i], 1, q1[i], 1 ); 
            
            Vmath::Svtvp(qBar1[i].num_elements(), F21, q0[i], 1, qBar1[i], 1, qBar1[i], 1 );    
            Vmath::Svtvp(qBar1[i].num_elements(), F22, qBar0[i], 1, qBar1[i], 1, qBar1[i], 1 );             
        }
        
        
        void DriverSteadyState::ConvergenceHistory(const Array<OneD, const Array<OneD, NekDouble> > &qBar1,
                                                   const Array<OneD, const Array<OneD, NekDouble> > &q0,
                                                   NekDouble &MaxNormDiff_q_qBar,
                                                   NekDouble &MaxNormDiff_q1_q0)
        {
            //This routine evaluates |q-qBar|_L2 and save the value in "ConvergenceHistory.txt"
Bastien Jordi's avatar
Bastien Jordi committed
232
            
Daniele de Grazia's avatar
Daniele de Grazia committed
233
234
            Array<OneD, NekDouble > NormDiff_q_qBar(NumVar_SFD, 1.0);
            Array<OneD, NekDouble > NormDiff_q1_q0(NumVar_SFD, 1.0);
Bastien Jordi's avatar
Bastien Jordi committed
235
236
237
238
239
            
            MaxNormDiff_q_qBar=0.0;
            MaxNormDiff_q1_q0=0.0;
            
            //Norm Calculation
Daniele de Grazia's avatar
Daniele de Grazia committed
240
            for(int i = 0; i < NumVar_SFD; ++i)
Bastien Jordi's avatar
Bastien Jordi committed
241
            {
Bastien Jordi's avatar
Bastien Jordi committed
242
243
244
                //To check convergence of SFD
                //NormDiff_q_qBar[i] = m_equ[0]->L2Error(i, qBar1[i], false);
                NormDiff_q_qBar[i] = m_equ[0]->LinfError(i, qBar1[i]);
Bastien Jordi's avatar
Bastien Jordi committed
245
                
Bastien Jordi's avatar
Bastien Jordi committed
246
247
                //To check convergence of Navier-Stokes
                NormDiff_q1_q0[i] = m_equ[0]->LinfError(i, q0[i]);
Bastien Jordi's avatar
Bastien Jordi committed
248
                
Bastien Jordi's avatar
Bastien Jordi committed
249
                if (MaxNormDiff_q_qBar < NormDiff_q_qBar[i])
250
                {
Bastien Jordi's avatar
Bastien Jordi committed
251
                    MaxNormDiff_q_qBar = NormDiff_q_qBar[i];
252
                }
Bastien Jordi's avatar
Bastien Jordi committed
253
                
Bastien Jordi's avatar
Bastien Jordi committed
254
                if (MaxNormDiff_q1_q0 < NormDiff_q1_q0[i])
Bastien Jordi's avatar
Bastien Jordi committed
255
                {
Bastien Jordi's avatar
Bastien Jordi committed
256
                    MaxNormDiff_q1_q0 = NormDiff_q1_q0[i];
Bastien Jordi's avatar
Bastien Jordi committed
257
                }
Bastien Jordi's avatar
Bastien Jordi committed
258
259
260
261
262
263
264
            }
            
            #if NEKTAR_USE_MPI   
            MPI_Comm_rank(MPI_COMM_WORLD,&MPIrank);
            if (MPIrank==0)
            {
                cout << "SFD (MPI) - Step: " << m_n+1 << "; Time: " << m_equ[0]->GetFinalTime() <<  "; |q-qBar|inf = " << MaxNormDiff_q_qBar << "; |q1-q0|inf = " << MaxNormDiff_q1_q0 << ";\t for X = " << m_X0 <<" and Delta = " << m_Delta0 <<endl;
265
                std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app); 
266
                m_file << m_n+1 << "\t" << m_equ[0]->GetFinalTime() << "\t" << MaxNormDiff_q_qBar << "\t" << MaxNormDiff_q1_q0 << endl;
267
268
                m_file.close();
            }
Bastien Jordi's avatar
Bastien Jordi committed
269
270
271
272
273
274
275
            #else
            cout << "SFD - Step: " << m_n+1 << "; Time: " << m_equ[0]->GetFinalTime() <<  "; |q-qBar|inf = " << MaxNormDiff_q_qBar << "; |q1-q0|inf = " << MaxNormDiff_q1_q0 << ";\t for X = " << m_X0 <<" and Delta = " << m_Delta0 <<endl;
            std::ofstream m_file( "ConvergenceHistory.txt", std::ios_base::app); 
            m_file << m_n+1 << "\t" << m_equ[0]->GetFinalTime() << "\t" << MaxNormDiff_q_qBar << "\t" << MaxNormDiff_q1_q0 << endl;
            m_file.close();
            #endif              
            
Bastien Jordi's avatar
Bastien Jordi committed
276
277
        }
    }
Bastien Jordi's avatar
Bastien Jordi committed
278
279
}