IncNavierStokes.cpp 21.1 KB
Newer Older
Spencer Sherwin's avatar
Spencer Sherwin committed
1 2
///////////////////////////////////////////////////////////////////////////////
//
Spencer Sherwin's avatar
Spencer Sherwin committed
3
// File IncNavierStokes.cpp
Spencer Sherwin's avatar
Spencer Sherwin committed
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
//
// 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 class definition built on
// ADRBase class
//
///////////////////////////////////////////////////////////////////////////////

37 38 39
#include <iomanip>
#include <boost/algorithm/string.hpp>

40
#include <IncNavierStokesSolver/EquationSystems/IncNavierStokes.h>
Spencer Sherwin's avatar
Spencer Sherwin committed
41
#include <LibUtilities/BasicUtils/Timer.h>
42 43
#include <LibUtilities/Communication/Comm.h>
#include <SolverUtils/Filters/Filter.h>
44 45
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/Expansion3D.h>
Spencer Sherwin's avatar
Spencer Sherwin committed
46 47 48

namespace Nektar
{
49

Spencer Sherwin's avatar
Spencer Sherwin committed
50 51 52 53 54 55
    /**
     * Constructor. Creates ...
     *
     * \param 
     * \param
     */
56
    IncNavierStokes::IncNavierStokes(const LibUtilities::SessionReaderSharedPtr& pSession):
57
        UnsteadySystem(pSession),
58
        m_subSteppingScheme(false),
Dave Moxey's avatar
Dave Moxey committed
59 60
        m_SmoothAdvection(false),
        m_steadyStateSteps(0)
Spencer Sherwin's avatar
Spencer Sherwin committed
61
    {
62 63 64 65
    }

    void IncNavierStokes::v_InitObject()
    {
66
        UnsteadySystem::v_InitObject();
67
        
Chris Cantwell's avatar
Chris Cantwell committed
68
        int i,j;
69
        int numfields = m_fields.num_elements();
70
        std::string velids[] = {"u","v","w"};
Spencer Sherwin's avatar
Spencer Sherwin committed
71
        
72 73
        // Set up Velocity field to point to the first m_expdim of m_fields; 
        m_velocity = Array<OneD,int>(m_spacedim);
Spencer Sherwin's avatar
Spencer Sherwin committed
74
        
75
        for(i = 0; i < m_spacedim; ++i)
Spencer Sherwin's avatar
Spencer Sherwin committed
76
        {
77
            for(j = 0; j < numfields; ++j)
Spencer Sherwin's avatar
Spencer Sherwin committed
78
            {
79
                std::string var = m_boundaryConditions->GetVariable(j);
80
                if(boost::iequals(velids[i], var))
81 82 83 84 85
                {
                    m_velocity[i] = j;
                    break;
                }
                
Chris Cantwell's avatar
Chris Cantwell committed
86
                ASSERTL0(j != numfields, "Failed to find field: " + var);
Spencer Sherwin's avatar
Spencer Sherwin committed
87 88
            }
        }
89
        
Spencer Sherwin's avatar
Spencer Sherwin committed
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
        // Set up equation type enum using kEquationTypeStr
        for(i = 0; i < (int) eEquationTypeSize; ++i)
        {
            bool match;
            m_session->MatchSolverInfo("EQTYPE",kEquationTypeStr[i],match,false);
            if(match)
            {
                m_equationType = (EquationType)i; 
                break;
            }
        }
        ASSERTL0(i != eEquationTypeSize,"EQTYPE not found in SOLVERINFO section");
        
        // This probably should to into specific implementations 
        // Equation specific Setups 
        switch(m_equationType)
        {
        case eSteadyStokes: 
        case eSteadyOseen: 
109
        case eSteadyNavierStokes:
110
        case eSteadyLinearisedNS: 
Spencer Sherwin's avatar
Spencer Sherwin committed
111 112 113
            break;
        case eUnsteadyNavierStokes:
        case eUnsteadyStokes:
114 115
            {
                m_session->LoadParameter("IO_InfoSteps", m_infosteps, 0);
116
                m_session->LoadParameter("IO_CFLSteps", m_cflsteps, 0);
117 118
                m_session->LoadParameter("SteadyStateSteps", m_steadyStateSteps, 0);
                m_session->LoadParameter("SteadyStateTol", m_steadyStateTol, 1e-6);
Spencer Sherwin's avatar
Spencer Sherwin committed
119
            
Gabriele Rocco's avatar
Gabriele Rocco committed
120
				
121 122 123 124
                // check to see if any user defined boundary condition is
                // indeed implemented
                
                for(int n = 0; n < m_fields[0]->GetBndConditions().num_elements(); ++n)
Chris Cantwell's avatar
Chris Cantwell committed
125
                {    
126 127
                    // Time Dependent Boundary Condition (if no user
                    // defined then this is empty)
Chris Cantwell's avatar
Chris Cantwell committed
128 129 130 131 132 133 134 135 136 137 138 139
                    ASSERTL0 (
                        m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
                            SpatialDomains::eNoUserDefined ||
                        m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
                            SpatialDomains::eWall_Forces ||
                        m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
                            SpatialDomains::eTimeDependent ||
                        m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
                            SpatialDomains::eRadiation ||
                        m_fields[0]->GetBndConditions()[n]->GetUserDefined() ==
                            SpatialDomains::eI,
                        "Unknown USERDEFINEDTYPE boundary condition");
Spencer Sherwin's avatar
Spencer Sherwin committed
140 141 142 143 144 145 146 147 148 149
                }
            }
            break;
        case eNoEquationType:
        default:
            ASSERTL0(false,"Unknown or undefined equation type");
        }
        
        m_session->LoadParameter("Kinvis", m_kinvis);
        
Chris Cantwell's avatar
Chris Cantwell committed
150 151 152
        // Default advection type per solver
        std::string vConvectiveType;
        switch(m_equationType)
153
        {
Chris Cantwell's avatar
Chris Cantwell committed
154 155 156 157 158 159 160 161 162 163
            case eUnsteadyStokes:
                vConvectiveType = "NoAdvection";
                break;
            case eUnsteadyNavierStokes:
            case eSteadyNavierStokes:
                vConvectiveType = "Convective";
                break;
            case eUnsteadyLinearisedNS:
                vConvectiveType = "Linearised";
                break;
164 165
            default:
                break;
166
        }
Chris Cantwell's avatar
Chris Cantwell committed
167 168 169

        // Check if advection type overridden
        if (m_session->DefinesTag("AdvectiveType") && m_equationType != eUnsteadyStokes)
170
        {
Chris Cantwell's avatar
Chris Cantwell committed
171
            vConvectiveType = m_session->GetTag("AdvectiveType");
172
        }
Chris Cantwell's avatar
Chris Cantwell committed
173 174 175

        // Initialise advection
        m_advObject = GetAdvectionTermFactory().CreateInstance(vConvectiveType, m_session, m_graph);
176
        
177
        // Forcing terms
178 179
        m_forcing = SolverUtils::Forcing::Load(m_session, m_fields,
                                               v_GetForceDimension());
180

181 182
        // check to see if any Robin boundary conditions and if so set
        // up m_field to boundary condition maps;
Chris Cantwell's avatar
Chris Cantwell committed
183 184 185
        m_fieldsBCToElmtID  = Array<OneD, Array<OneD, int> >(numfields);
        m_fieldsBCToTraceID = Array<OneD, Array<OneD, int> >(numfields);
        m_fieldsRadiationFactor  = Array<OneD, Array<OneD, NekDouble> > (numfields);
186 187
        
        for (i = 0; i < m_fields.num_elements(); ++i)
188
        {
189 190 191 192 193 194 195 196 197
            bool Set = false;

            Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
            Array<OneD, MultiRegions::ExpListSharedPtr>  BndExp;
            int radpts = 0;
            
            BndConds = m_fields[i]->GetBndConditions();
            BndExp   = m_fields[i]->GetBndCondExpansions();
            for(int n = 0; n < BndConds.num_elements(); ++n)
Chris Cantwell's avatar
Chris Cantwell committed
198
            {    
199 200
                if(BndConds[n]->GetUserDefined() == SpatialDomains::eRadiation)
                {
201 202 203
                    ASSERTL0(BndConds[n]->GetBoundaryConditionType() == SpatialDomains::eRobin,
                             "Radiation boundary condition must be of type Robin <R>");
                    
204 205 206 207 208 209 210 211 212 213 214 215 216 217
                    if(Set == false)
                    {
                        m_fields[i]->GetBoundaryToElmtMap(m_fieldsBCToElmtID[i],m_fieldsBCToTraceID[i]);
                        Set = true;
                    }
                    radpts += BndExp[n]->GetTotPoints();
                }
            }

            m_fieldsRadiationFactor[i] = Array<OneD, NekDouble>(radpts);

            radpts = 0; // reset to use as a counter

            for(int n = 0; n < BndConds.num_elements(); ++n)
Chris Cantwell's avatar
Chris Cantwell committed
218
            {    
219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
                if(BndConds[n]->GetUserDefined() == SpatialDomains::eRadiation)
                {
                    
                    int npoints    = BndExp[n]->GetNpoints();
                    Array<OneD, NekDouble> x0(npoints,0.0);
                    Array<OneD, NekDouble> x1(npoints,0.0);
                    Array<OneD, NekDouble> x2(npoints,0.0);
                    Array<OneD, NekDouble> tmpArray;

                    BndExp[n]->GetCoords(x0,x1,x2);
                    
                    LibUtilities::Equation coeff = 
                        boost::static_pointer_cast<
                    SpatialDomains::RobinBoundaryCondition
                        >(BndConds[n])->m_robinPrimitiveCoeff;
                    
                    coeff.Evaluate(x0,x1,x2,m_time, 
                                   tmpArray = m_fieldsRadiationFactor[i]+ radpts);
                    //Vmath::Neg(npoints,tmpArray = m_fieldsRadiationFactor[i]+ radpts,1);
                    radpts += npoints;
                }
            }
241
        }
242 243

        // Set up Field Meta Data for output files
244 245
        m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
        m_fieldMetaDataMap["TimeStep"] = boost::lexical_cast<std::string>(m_timestep);
246
		
Andrew Comerford's avatar
Andrew Comerford committed
247 248 249
        // creation of the extrapolation object
        if(m_equationType == eUnsteadyNavierStokes)
        {
250 251 252
            std::string vExtrapolation = "Standard";

            if (m_session->DefinesSolverInfo("Extrapolation"))
253
            {
254
                vExtrapolation = m_session->GetSolverInfo("Extrapolation");
255
            }
Andrew Comerford's avatar
Andrew Comerford committed
256 257 258 259 260 261 262 263
                        
            m_extrapolation = GetExtrapolateFactory().CreateInstance(
                vExtrapolation, 
                m_session,
                m_fields,
                m_velocity,
                m_advObject);
        }
264 265
    }

Andrew Comerford's avatar
Andrew Comerford committed
266 267 268
    /**
     * Destructor
     */
269
    IncNavierStokes::~IncNavierStokes(void)
270
    {
Spencer Sherwin's avatar
Spencer Sherwin committed
271
    }
272 273

    
274
	/**
Ale's avatar
Ale committed
275
	 *
276
	 */
277 278 279 280 281 282 283 284 285 286 287 288
    void IncNavierStokes::v_GetFluxVector(const int i, 
                                          Array<OneD, Array<OneD, NekDouble> > &physfield,
                                            Array<OneD, Array<OneD, NekDouble> > &flux)
    {
        ASSERTL1(flux.num_elements() == m_velocity.num_elements(),"Dimension of flux array and velocity array do not match");

        for(int j = 0; j < flux.num_elements(); ++j)
        {
            Vmath::Vmul(GetNpoints(), physfield[i], 1, m_fields[m_velocity[j]]->GetPhys(), 1, flux[j], 1);
        }
    }

289 290 291
    /**
	 * Calcualate numerical fluxes
	 */
292 293 294 295 296 297 298 299 300 301 302 303 304
    void IncNavierStokes::v_NumericalFlux(Array<OneD, Array<OneD, NekDouble> > &physfield, 
                                          Array<OneD, Array<OneD, NekDouble> > &numflux)
    {
        /// Counter variable
        int i;

        /// Number of trace points
        int nTracePts   = GetTraceNpoints();
        
        /// Number of spatial dimensions
        int nDimensions = m_spacedim;

        /// Forward state array
305
        Array<OneD, NekDouble> Fwd(2*nTracePts);
306 307
        
        /// Backward state array
308
        Array<OneD, NekDouble> Bwd = Fwd + nTracePts;
309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
        
        /// Normal velocity array
        Array<OneD, NekDouble> Vn (nTracePts, 0.0);
        
        // Extract velocity field along the trace space and multiply by trace normals
        for(i = 0; i < nDimensions; ++i)
        {
            m_fields[0]->ExtractTracePhys(m_fields[m_velocity[i]]->GetPhys(), Fwd);
            Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, Fwd, 1, Vn, 1, Vn, 1);
        }

        /// Compute the numerical fluxes at the trace points
        for(i = 0; i < numflux.num_elements(); ++i)
        {
            /// Extract forwards/backwards trace spaces
            m_fields[i]->GetFwdBwdTracePhys(physfield[i], Fwd, Bwd);

            /// Upwind between elements
            m_fields[i]->GetTrace()->Upwind(Vn, Fwd, Bwd, numflux[i]);

329 330
            /// Calculate the numerical fluxes multipling Fwd or Bwd
            /// by the normal advection velocity
331 332 333
            Vmath::Vmul(nTracePts, numflux[i], 1, Vn, 1, numflux[i], 1);
        }
    }
334

335 336 337
	/**
	 * Evaluation -N(V) for all fields except pressure using m_velocity
	 */
338 339 340
    void IncNavierStokes::EvaluateAdvectionTerms(const Array<OneD, const Array<OneD, NekDouble> > &inarray, 
                                                 Array<OneD, Array<OneD, NekDouble> > &outarray, 
                                                 Array<OneD, NekDouble> &wk)
Spencer Sherwin's avatar
Spencer Sherwin committed
341
    {
Chris Cantwell's avatar
Chris Cantwell committed
342
        int i;
343 344 345 346 347 348 349 350
        int nqtot      = m_fields[0]->GetTotPoints();
        int VelDim     = m_velocity.num_elements();
        Array<OneD, Array<OneD, NekDouble> > velocity(VelDim);
        Array<OneD, NekDouble > Deriv;
        for(i = 0; i < VelDim; ++i)
        {
            velocity[i] = inarray[m_velocity[i]]; 
        }
Spencer Sherwin's avatar
Spencer Sherwin committed
351 352 353
        // Set up Derivative work space; 
        if(wk.num_elements())
        {
Chris Cantwell's avatar
Chris Cantwell committed
354 355
            ASSERTL0(wk.num_elements() >= nqtot*VelDim,
                     "Workspace is not sufficient");
Spencer Sherwin's avatar
Spencer Sherwin committed
356 357 358 359
            Deriv = wk;
        }
        else
        {
360
            Deriv = Array<OneD, NekDouble> (nqtot*VelDim);
Spencer Sherwin's avatar
Spencer Sherwin committed
361
        }
362

363
        m_advObject->DoAdvection(m_fields,m_nConvectiveFields, 
364
                                 m_velocity,inarray,outarray,m_time,Deriv);
Spencer Sherwin's avatar
Spencer Sherwin committed
365
    }
366
    
367 368 369
    /**
	 * Time dependent boundary conditions updating
	 */
370 371
    void IncNavierStokes::SetBoundaryConditions(NekDouble time)
    {
Chris Cantwell's avatar
Chris Cantwell committed
372
        int i, n;
373 374
        std::string varName;
        int nvariables = m_fields.num_elements();
375
        
Chris Cantwell's avatar
Chris Cantwell committed
376
        for (i = 0; i < nvariables; ++i)
377
        {
Chris Cantwell's avatar
Chris Cantwell committed
378
            for(n = 0; n < m_fields[i]->GetBndConditions().num_elements(); ++n)
Chris Cantwell's avatar
Chris Cantwell committed
379
            {    
380 381
                if(m_fields[i]->GetBndConditions()[n]->GetUserDefined() ==
                   SpatialDomains::eTimeDependent)
382
                {
383
                    varName = m_session->GetVariable(i);
384
                    m_fields[i]->EvaluateBoundaryConditions(time, varName);
385
                }
386

387
            }
388

389
            // Set Radiation conditions if required
Chris Cantwell's avatar
Chris Cantwell committed
390
            SetRadiationBoundaryForcing(i);
391 392 393
        }
    }
    
394 395 396
    /**
	 * Probably should be pushed back into ContField? 
	 */
397 398 399 400 401 402 403 404 405 406 407 408
    void IncNavierStokes::SetRadiationBoundaryForcing(int fieldid)
    {
        int  i,n;
        
        Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
        Array<OneD, MultiRegions::ExpListSharedPtr>  BndExp;
        
        
        BndConds = m_fields[fieldid]->GetBndConditions();
        BndExp   = m_fields[fieldid]->GetBndCondExpansions();
        
        StdRegions::StdExpansionSharedPtr   elmt;
409
        StdRegions::StdExpansionSharedPtr Bc;
410 411 412 413 414 415
        
        int cnt;
        int elmtid,nq,offset, boundary;
        Array<OneD, NekDouble> Bvals, U;
        int cnt1 = 0;

416

417 418 419 420 421 422 423 424 425 426 427 428 429
        for(cnt = n = 0; n < BndConds.num_elements(); ++n)
        {            
            SpatialDomains::BndUserDefinedType type = BndConds[n]->GetUserDefined(); 
            
            if((BndConds[n]->GetBoundaryConditionType() == SpatialDomains::eRobin)&&(type == SpatialDomains::eRadiation))
            {
                for(i = 0; i < BndExp[n]->GetExpSize(); ++i,cnt++)
                {
                    elmtid = m_fieldsBCToElmtID[fieldid][cnt];
                    elmt   = m_fields[fieldid]->GetExp(elmtid);
                    offset = m_fields[fieldid]->GetPhys_Offset(elmtid);
                    
                    U = m_fields[fieldid]->UpdatePhys() + offset;
430
                    Bc = BndExp[n]->GetExp(i);
431 432 433 434 435 436
                    
                    boundary = m_fieldsBCToTraceID[fieldid][cnt];
                    
                    // Get edge values and put into ubc
                    nq = Bc->GetTotPoints();
                    Array<OneD, NekDouble> ubc(nq);
437
                    elmt->GetTracePhysVals(boundary,Bc,U,ubc);
438
                    
439 440
                    Vmath::Vmul(nq,&m_fieldsRadiationFactor[fieldid][cnt1 + 
                             BndExp[n]->GetPhys_Offset(i)],1,&ubc[0],1,&ubc[0],1);
441 442 443 444 445 446 447

                    Bvals = BndExp[n]->UpdateCoeffs()+BndExp[n]->GetCoeff_Offset(i);

                    Bc->IProductWRTBase(ubc,Bvals); 
                }
                cnt1 += BndExp[n]->GetTotPoints();
            }
Chris Cantwell's avatar
Chris Cantwell committed
448
            else if(type == SpatialDomains::eNoUserDefined || type == SpatialDomains::eWall_Forces || type == SpatialDomains::eTimeDependent || type == SpatialDomains::eHigh) 
449 450 451 452 453 454 455 456 457 458
            {
                cnt += BndExp[n]->GetExpSize();
            }
            else
            {
                ASSERTL0(false,"Unknown USERDEFINEDTYPE in pressure boundary condition");
            }
        }
    }

459

Chris Cantwell's avatar
Chris Cantwell committed
460 461 462
    /**
     * Add an additional forcing term programmatically.
     */
463 464 465 466 467 468
    void IncNavierStokes::AddForcing(const SolverUtils::ForcingSharedPtr& pForce)
    {
        m_forcing.push_back(pForce);
    }


Chris Cantwell's avatar
Chris Cantwell committed
469 470
    /**
     * Decide if at a steady state if the discrerte L2 sum of the
471 472 473
	 * coefficients is the same as the previous step to within the
	 * tolerance m_steadyStateTol;
	 */
474 475 476 477
    bool IncNavierStokes::CalcSteadyState(void)
    {
        static NekDouble previousL2 = 0.0;
        bool returnval = false;
478
        
479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497
        NekDouble L2 = 0.0;
        
        // calculate L2 discrete summation 
        int ncoeffs = m_fields[0]->GetNcoeffs(); 
        
        for(int i = 0; i < m_fields.num_elements(); ++i)
        {
            L2 += Vmath::Dot(ncoeffs,m_fields[i]->GetCoeffs(),1,m_fields[i]->GetCoeffs(),1);
        }
        
        if(fabs(L2-previousL2) < ncoeffs*m_steadyStateTol)
        {
            returnval = true;
        }

        previousL2 = L2;

        return returnval;
    }
Gianmarco Mengaldo's avatar
Gianmarco Mengaldo committed
498
    
499 500 501
	/**
	 *
	 */
502
    Array<OneD, NekDouble> IncNavierStokes::GetElmtCFLVals(void)
Chris Cantwell's avatar
Chris Cantwell committed
503 504
    {
        int n_vel     = m_velocity.num_elements();
Dave Moxey's avatar
Dave Moxey committed
505
        int n_element = m_fields[0]->GetExpSize(); 
506 507 508 509
        
        const Array<OneD, int> ExpOrder = GetNumExpModesPerExp();
        Array<OneD, int> ExpOrderList (n_element, ExpOrder);
        
Dave Moxey's avatar
Dave Moxey committed
510
        const NekDouble cLambda = 0.2; // Spencer book pag. 317
511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526
        
        Array<OneD, NekDouble> cfl        (n_element, 0.0);
        Array<OneD, NekDouble> stdVelocity(n_element, 0.0);
        Array<OneD, Array<OneD, NekDouble> > velfields; 
        
        if(m_HomogeneousType == eHomogeneous1D) // just do check on 2D info
        {
            velfields = Array<OneD, Array<OneD, NekDouble> >(2);

            for(int i = 0; i < 2; ++i)
            {
                velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
            }        
        }
        else
        {
Chris Cantwell's avatar
Chris Cantwell committed
527
            velfields = Array<OneD, Array<OneD, NekDouble> >(n_vel);
528

Chris Cantwell's avatar
Chris Cantwell committed
529
            for(int i = 0; i < n_vel; ++i)
530 531 532 533
            {
                velfields[i] = m_fields[m_velocity[i]]->UpdatePhys();
            }        
        }
534 535

        stdVelocity = m_extrapolation->GetMaxStdVelocity(velfields);
536 537 538 539 540 541 542
        
        for(int el = 0; el < n_element; ++el)
        {
            cfl[el] =  m_timestep*(stdVelocity[el] * cLambda *
                                   (ExpOrder[el]-1) * (ExpOrder[el]-1));
        }
        
543 544 545
        return cfl;
    }
    
546 547 548
    /**
	 *
	 */
549 550 551 552 553 554
    NekDouble IncNavierStokes::GetCFLEstimate(int &elmtid)
    { 
        int n_element = m_fields[0]->GetExpSize(); 

        Array<OneD, NekDouble> cfl = GetElmtCFLVals();
        
555
        elmtid = Vmath::Imax(n_element,cfl,1);
556 557 558 559
        NekDouble CFL,CFL_loc;

        CFL = CFL_loc = cfl[elmtid];
        m_comm->AllReduce(CFL,LibUtilities::ReduceMax);
560 561

        // unshuffle elmt id if data is not stored in consecutive order. 
562
        elmtid = m_fields[0]->GetExp(elmtid)->GetGeom()->GetGlobalID();
563 564 565 566
        if(CFL != CFL_loc)
        {
            elmtid = -1;
        }
567

568
        m_comm->AllReduce(elmtid,LibUtilities::ReduceMax);
569
        
Chris Cantwell's avatar
Chris Cantwell committed
570 571
        // express element id with respect to plane
        if(m_HomogeneousType == eHomogeneous1D)
572 573 574
        {
            elmtid = elmtid%m_fields[0]->GetPlane(0)->GetExpSize();
        }
575 576
        return CFL;
    }
577

Chris Cantwell's avatar
Chris Cantwell committed
578 579 580 581

    /**
     * Perform the extrapolation.
     */
582 583 584 585 586 587 588
    bool IncNavierStokes::v_PreIntegrate(int step)
    {
        m_extrapolation->SubStepSaveFields(step);
        m_extrapolation->SubStepAdvance(m_intSoln,step,m_time);
        return false;
    }

Chris Cantwell's avatar
Chris Cantwell committed
589 590 591 592

    /**
     * Estimate CFL and perform steady-state check
     */
593 594 595 596 597 598 599 600 601
    bool IncNavierStokes::v_PostIntegrate(int step)
    {
        if(m_cflsteps && !((step+1)%m_cflsteps))
        {
            int elmtid;
            NekDouble cfl = GetCFLEstimate(elmtid);

            if(m_comm->GetRank() == 0)
            {
Chris Cantwell's avatar
Chris Cantwell committed
602 603
                cout << "CFL (zero plane): "<< cfl << " (in elmt "
                     << elmtid << ")" << endl;
604 605 606 607 608 609 610
            }
        }

        if(m_steadyStateSteps && step && (!((step+1)%m_steadyStateSteps)))
        {
            if(CalcSteadyState() == true)
            {
Chris Cantwell's avatar
Chris Cantwell committed
611 612
                cout << "Reached Steady State to tolerance "
                     << m_steadyStateTol << endl;
613 614 615 616 617 618
                return true;
            }
        }

        return false;
    }
Spencer Sherwin's avatar
Spencer Sherwin committed
619 620
} //end of namespace