VCSMapping.cpp 36.5 KB
Newer Older
1 2
///////////////////////////////////////////////////////////////////////////////
//
Douglas Serson's avatar
Douglas Serson committed
3
// File VCSMapping.cpp
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
//
// 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: Velocity Correction Scheme w/ coordinate transformation 
// for the Incompressible Navier Stokes equations
///////////////////////////////////////////////////////////////////////////////

#include <IncNavierStokesSolver/EquationSystems/VCSMapping.h>
#include <LibUtilities/BasicUtils/Timer.h>
#include <SolverUtils/Core/Misc.h>

#include <boost/algorithm/string.hpp>

Chris Cantwell's avatar
Chris Cantwell committed
42 43
using namespace std;

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
namespace Nektar
{
    string VCSMapping::className = 
        SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
            "VCSMapping", 
            VCSMapping::create);

    /**
     * Constructor. Creates ...
     *
     * \param 
     * \param
     */
    VCSMapping::VCSMapping(
            const LibUtilities::SessionReaderSharedPtr& pSession)
        : UnsteadySystem(pSession),
          VelocityCorrectionScheme(pSession)  
    {

    }

    void VCSMapping::v_InitObject()
    {
        VelocityCorrectionScheme::v_InitObject();
        
69
        m_mapping = GlobalMapping::Mapping::Load(m_session, m_fields); 
70
        ASSERTL0(m_mapping,
Douglas Serson's avatar
Douglas Serson committed
71
             "Could not create mapping in VCSMapping.");
72 73 74 75 76 77 78 79 80
        
        std::string vExtrapolation = "Mapping";
        m_extrapolation = GetExtrapolateFactory().CreateInstance(
            vExtrapolation,
            m_session,
            m_fields,
            m_pressure,
            m_velocity,
            m_advObject); 
Douglas Serson's avatar
Douglas Serson committed
81 82
        m_extrapolation->SubSteppingTimeIntegration(
                            m_intScheme->GetIntegrationMethod(), m_intScheme);
83
        m_extrapolation->GenerateHOPBCMap();        
84 85 86

       // Storage to extrapolate pressure forcing
        int physTot = m_fields[0]->GetTotPoints();
Douglas Serson's avatar
Douglas Serson committed
87
        int intSteps = 1;
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
        int intMethod = m_intScheme->GetIntegrationMethod();
        switch(intMethod)
        {
            case LibUtilities::eIMEXOrder1:
            {
                intSteps = 1; 
            }
            break;
            case LibUtilities::eIMEXOrder2:
            {
                intSteps = 2;
            }
            break;
            case LibUtilities::eIMEXOrder3:
            {
                intSteps = 3;
            }
            break;
        }        
Douglas Serson's avatar
Douglas Serson committed
107
        m_presForcingCorrection= Array<OneD, Array<OneD, NekDouble> >(intSteps);
108 109 110 111
        for(int i = 0; i < m_presForcingCorrection.num_elements(); i++)
        {
            m_presForcingCorrection[i] = Array<OneD, NekDouble>(physTot,0.0);
        }
112
        m_verbose = (m_session->DefinesCmdLineArgument("verbose"))? true :false;
113 114 115
        
        // Load solve parameters related to the mapping
        // Flags determining if pressure/viscous terms should be treated implicitly
Douglas Serson's avatar
Douglas Serson committed
116 117 118 119 120 121
        m_session->MatchSolverInfo("MappingImplicitPressure","True",
                                                    m_implicitPressure,false);
        m_session->MatchSolverInfo("MappingImplicitViscous","True",
                                                    m_implicitViscous,false);
        m_session->MatchSolverInfo("MappingNeglectViscous","True",
                                                    m_neglectViscous,false);
122 123 124 125 126
        
        if (m_neglectViscous)
        {
            m_implicitViscous = false;
        }
127 128
        
        // Tolerances and relaxation parameters for implicit terms
Douglas Serson's avatar
Douglas Serson committed
129 130 131 132 133 134 135 136
        m_session->LoadParameter("MappingPressureTolerance",
                                            m_pressureTolerance,1e-12);
        m_session->LoadParameter("MappingViscousTolerance",
                                            m_viscousTolerance,1e-12);
        m_session->LoadParameter("MappingPressureRelaxation",
                                            m_pressureRelaxation,1.0);
        m_session->LoadParameter("MappingViscousRelaxation",
                                            m_viscousRelaxation,1.0);
137
        
138 139 140 141 142 143 144 145 146
    }
    
    /**
     * Destructor
     */
    VCSMapping::~VCSMapping(void)
    {        
    }
    
147 148 149
    void VCSMapping::v_DoInitialise(void)
    {
        UnsteadySystem::v_DoInitialise();
150

151
        // Set up Field Meta Data for output files
Douglas Serson's avatar
Douglas Serson committed
152 153 154 155
        m_fieldMetaDataMap["Kinvis"]   = 
                                boost::lexical_cast<std::string>(m_kinvis);
        m_fieldMetaDataMap["TimeStep"] = 
                                boost::lexical_cast<std::string>(m_timestep);
156 157

        // Correct Dirichlet boundary conditions to account for mapping
158
        m_mapping->UpdateBCs(0.0);
159 160 161 162 163 164 165 166 167
        //
        for(int i = 0; i < m_nConvectiveFields; ++i)
        {
            m_fields[i]->LocalToGlobal();
            m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
            m_fields[i]->GlobalToLocal();
            m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
                                  m_fields[i]->UpdatePhys());
        }
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182

        // Initialise m_gradP
        int physTot = m_fields[0]->GetTotPoints();
        m_gradP = Array<OneD, Array<OneD, NekDouble> >(m_nConvectiveFields);
        for(int i = 0; i < m_nConvectiveFields; ++i)
        {
            m_gradP[i] = Array<OneD, NekDouble>(physTot,0.0);
            m_pressure->PhysDeriv(MultiRegions::DirCartesianMap[i],
                                    m_pressure->GetPhys(),
                                    m_gradP[i]);
            if(m_pressure->GetWaveSpace())
            {
                m_pressure->HomogeneousBwdTrans(m_gradP[i],m_gradP[i]);
            }
        }
183
    }
184 185 186 187 188 189 190 191

    /**
     * Explicit part of the method - Advection, Forcing + HOPBCs
     */
    void VCSMapping::v_EvaluateAdvection_SetPressureBCs(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray, 
        Array<OneD, Array<OneD, NekDouble> > &outarray, 
        const NekDouble time)
192
    {       
193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
        EvaluateAdvectionTerms(inarray, outarray);

        // Smooth advection
        if(m_SmoothAdvection)
        {
            for(int i = 0; i < m_nConvectiveFields; ++i)
            {
                m_pressure->SmoothField(outarray[i]);
            }
        }

        // Add forcing terms
        std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
        for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
        {
            (*x)->Apply(m_fields, inarray, outarray, time);
        }
        
        // Add mapping terms
212
        ApplyIncNSMappingForcing( inarray, outarray);
213

214 215
        // Calculate High-Order pressure boundary conditions
        m_extrapolation->EvaluatePressureBCs(inarray,outarray,m_kinvis);
216 217 218 219 220 221 222 223 224 225 226 227

        // Update mapping and deal with Dirichlet boundary conditions
        if (m_mapping->IsTimeDependent())
        {
            if (m_mapping->IsFromFunction())
            {
                // If the transformation is explicitly defined, update it here
                // Otherwise, it will be done somewhere else (ForcingMovingBody)
                m_mapping->UpdateMapping(time+m_timestep);
            }
            m_mapping->UpdateBCs(time+m_timestep);
        }
228 229 230 231 232 233 234 235 236 237 238
    }
    
        
    /**
     * Forcing term for Poisson solver solver
     */ 
    void   VCSMapping::v_SetUpPressureForcing(
        const Array<OneD, const Array<OneD, NekDouble> > &fields, 
        Array<OneD, Array<OneD, NekDouble> > &Forcing, 
        const NekDouble aii_Dt)
    {                
239
        if (m_mapping->HasConstantJacobian())
240
        {
Douglas Serson's avatar
Douglas Serson committed
241 242
            VelocityCorrectionScheme::v_SetUpPressureForcing(fields, 
                                                            Forcing, aii_Dt);
243
        }
244 245 246 247 248 249 250 251 252
        else
        {
            int physTot = m_fields[0]->GetTotPoints();
            int nvel = m_nConvectiveFields;
            Array<OneD, NekDouble> wk(physTot, 0.0);

            Array<OneD, NekDouble> Jac(physTot,0.0);
            m_mapping->GetJacobian(Jac);
            
253
            // Calculate div(J*u/Dt)
254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277
            Vmath::Zero(physTot,Forcing[0],1);
            for(int i = 0; i < nvel; ++i)
            {
                if (m_fields[i]->GetWaveSpace())
                {
                    m_fields[i]->HomogeneousBwdTrans(fields[i],wk);
                }
                else
                {
                    Vmath::Vcopy(physTot, fields[i], 1, wk, 1);
                }
                Vmath::Vmul(physTot,wk,1,Jac,1,wk,1);
                if (m_fields[i]->GetWaveSpace())
                {
                    m_fields[i]->HomogeneousFwdTrans(wk,wk);
                }               
                m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[i],wk, wk);
                Vmath::Vadd(physTot,wk,1,Forcing[0],1,Forcing[0],1);
            }
            Vmath::Smul(physTot,1.0/aii_Dt,Forcing[0],1,Forcing[0],1); 

            //
            // If the mapping viscous terms are being treated explicitly 
            //        we need to apply a correction to the forcing
278
            if (!m_implicitViscous)
279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320
            {               
                bool wavespace = m_fields[0]->GetWaveSpace();
                m_fields[0]->SetWaveSpace(false);

                // 
                //  Part 1: div(J*grad(U/J . grad(J)))
                Array<OneD, Array<OneD, NekDouble> > tmp (nvel);
                Array<OneD, Array<OneD, NekDouble> > velocity (nvel);
                for(int i = 0; i < tmp.num_elements(); i++)
                {
                    tmp[i] = Array<OneD, NekDouble>(physTot,0.0);
                    velocity[i] = Array<OneD, NekDouble>(physTot,0.0);
                    if (wavespace)
                    {
                        m_fields[0]->HomogeneousBwdTrans(m_fields[i]->GetPhys(),
                                                         velocity[i]);
                    }
                    else
                    {
                        Vmath::Vcopy(physTot, m_fields[i]->GetPhys(), 1, 
                                                velocity[i], 1);
                    }
                }
                // Calculate wk = U.grad(J)
                m_mapping->DotGradJacobian(velocity, wk);
                // Calculate wk = (U.grad(J))/J
                Vmath::Vdiv(physTot, wk, 1, Jac, 1, wk, 1);
                // J*grad[(U.grad(J))/J]
                for(int i = 0; i < nvel; ++i)
                {
                    m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], 
                                wk, tmp[i]);
                    Vmath::Vmul(physTot, Jac, 1, tmp[i], 1, tmp[i], 1);
                }          
                // div(J*grad[(U.grad(J))/J])
                Vmath::Zero(physTot, wk, 1);
                for(int i = 0; i < nvel; ++i)
                {
                    m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], 
                                tmp[i], tmp[i]);
                    Vmath::Vadd(physTot, wk, 1, tmp[i], 1, wk, 1);
                }        
321

322
                // Part 2: grad(J) . curl(curl(U))
323
                m_mapping->CurlCurlField(velocity, tmp, m_implicitViscous);
Douglas Serson's avatar
Douglas Serson committed
324 325
                // dont need velocity any more, so reuse it
                m_mapping->DotGradJacobian(tmp, velocity[0]); 
326 327 328 329 330 331 332 333

                // Add two parts
                Vmath::Vadd(physTot, velocity[0], 1, wk, 1, wk, 1);

                // Multiply by kinvis
                Vmath::Smul(physTot, m_kinvis, wk, 1, wk, 1);

                // Extrapolate correction
Douglas Serson's avatar
Douglas Serson committed
334 335
                m_extrapolation->ExtrapolateArray(m_presForcingCorrection, 
                                                    wk, wk);
336 337 338 339 340 341 342 343 344 345 346 347

                // Put in wavespace
                if (wavespace)
                {
                    m_fields[0]->HomogeneousFwdTrans(wk,wk);
                }
                // Apply correction: Forcing = Forcing - correction
                Vmath::Vsub(physTot, Forcing[0], 1, wk, 1, Forcing[0], 1);

                m_fields[0]->SetWaveSpace(wavespace); 
            }
        }
348 349 350 351 352 353 354 355 356 357
    }
    
    /**
     * Forcing term for Helmholtz solver
     */
    void   VCSMapping::v_SetUpViscousForcing(
        const Array<OneD, const Array<OneD, NekDouble> > &inarray, 
        Array<OneD, Array<OneD, NekDouble> > &Forcing, 
        const NekDouble aii_Dt)
    {
358 359 360 361 362 363 364 365
        NekDouble aii_dtinv = 1.0/aii_Dt;
        int physTot = m_fields[0]->GetTotPoints();

        // Grad p
        m_pressure->BwdTrans(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());

        int nvel = m_velocity.num_elements();
        if(nvel == 2)
366
        {
367
            m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1]);
368 369 370
        }
        else
        {
371 372 373
            m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1],
                                  Forcing[2]);
        }
374

375 376 377 378
        // Copy grad p in physical space to m_gradP to reuse later
        if (m_pressure->GetWaveSpace())
        {
            for (int i=0; i<nvel; i++)
379
            {
380
                m_pressure->HomogeneousBwdTrans(Forcing[i],m_gradP[i]);
381
            }
382 383 384 385
        }
        else
        {
            for (int i=0; i<nvel; i++)
386
            {
387
                Vmath::Vcopy(physTot, Forcing[i], 1, m_gradP[i], 1);
388
            }
389 390 391 392
        }

        if ( (!m_mapping->HasConstantJacobian()) || m_implicitPressure)
        {
Douglas Serson's avatar
Douglas Serson committed
393
            // If pressure terms are treated explicitly, we need to divide by J
394
            //    if they are implicit, we need to calculate G(p)
395
            if (m_implicitPressure)
396 397
            {
                m_mapping->RaiseIndex(m_gradP, Forcing);
398 399 400 401 402 403 404
            }
            else
            {
                Array<OneD, NekDouble> Jac(physTot,0.0);
                m_mapping->GetJacobian(Jac);
                for (int i=0; i<nvel; i++)
                {
405 406
                    Vmath::Vdiv(physTot, m_gradP[i], 1, Jac, 1, Forcing[i], 1);
                }
407 408 409 410 411 412 413
            }
            // Transform back to wavespace
            if (m_pressure->GetWaveSpace())
            {
                for (int i=0; i<nvel; i++)
                {
                    m_pressure->HomogeneousFwdTrans(Forcing[i],Forcing[i]);
414 415 416 417 418 419 420 421 422 423
                }
            }
        }

        // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
        // need to be updated for the convected fields.
        for(int i = 0; i < nvel; ++i)
        {
            Blas::Daxpy(physTot,-aii_dtinv,inarray[i],1,Forcing[i],1);
            Blas::Dscal(physTot,1.0/m_kinvis,&(Forcing[i])[0],1);
424 425 426 427 428 429 430 431 432
        }
    }
    
    /**
     * Solve pressure system
     */
    void   VCSMapping::v_SolvePressure(
        const Array<OneD, NekDouble>  &Forcing)
    {
433
        if (!m_implicitPressure)
434 435 436 437 438 439 440
        {
            VelocityCorrectionScheme::v_SolvePressure(Forcing);
        }
        else
        {
            int physTot = m_fields[0]->GetTotPoints();
            int nvel = m_nConvectiveFields;
Douglas Serson's avatar
Douglas Serson committed
441 442 443 444
            bool converged = false;          // flag to mark if system converged
            int s = 0;                       // iteration counter
            NekDouble error;                 // L2 error at current iteration
            NekDouble forcing_L2 = 0.0;      // L2 norm of F
Douglas Serson's avatar
Douglas Serson committed
445 446 447
            
            int maxIter;
            m_session->LoadParameter("MappingMaxIter",maxIter,5000);
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471

            // rhs of the equation at current iteration
            Array< OneD, NekDouble> F_corrected(physTot, 0.0);
            // Pressure field at previous iteration
            Array<OneD, NekDouble>  previous_iter (physTot, 0.0);
            // Temporary variables
            Array<OneD, Array<OneD, NekDouble> > wk1(nvel);
            Array<OneD, Array<OneD, NekDouble> > wk2(nvel);
            Array<OneD, Array<OneD, NekDouble> > gradP(nvel);
            for(int i = 0; i < nvel; ++i)
            {
                wk1[i]   = Array<OneD, NekDouble> (physTot, 0.0);
                wk2[i]   = Array<OneD, NekDouble> (physTot, 0.0);
                gradP[i] = Array<OneD, NekDouble> (physTot, 0.0);
            }

            // Jacobian
            Array<OneD, NekDouble> Jac(physTot, 0.0);
            m_mapping->GetJacobian(Jac);

            // Factors for Laplacian system
            StdRegions::ConstFactorMap factors;
            factors[StdRegions::eFactorLambda] = 0.0;

Douglas Serson's avatar
Douglas Serson committed
472 473
            m_pressure->BwdTrans(m_pressure->GetCoeffs(),
                                    m_pressure->UpdatePhys());
474 475 476 477 478
            forcing_L2 = m_pressure->L2(Forcing, wk1[0]);
            while (!converged)
            {
                // Update iteration counter and set previous iteration field
                // (use previous timestep solution for first iteration)
Douglas Serson's avatar
Douglas Serson committed
479 480 481 482
                s++;         
                ASSERTL0(s < maxIter,
                         "VCSMapping exceeded maximum number of iterations.");
                
Douglas Serson's avatar
Douglas Serson committed
483 484
                Vmath::Vcopy(physTot, m_pressure->GetPhys(), 1, 
                                        previous_iter, 1);
485 486
                
                // Correct pressure bc to account for iteration
Douglas Serson's avatar
Douglas Serson committed
487
                m_extrapolation->CorrectPressureBCs(previous_iter);
488 489 490 491 492 493

                //
                // Calculate forcing term for this iteration
                //
                for(int i = 0; i < nvel; ++i)
                {
Douglas Serson's avatar
Douglas Serson committed
494 495
                    m_pressure->PhysDeriv(MultiRegions::DirCartesianMap[i], 
                                                    previous_iter, gradP[i]);
496 497 498 499 500 501 502 503 504 505 506 507 508 509
                    if(m_pressure->GetWaveSpace())
                    {
                        m_pressure->HomogeneousBwdTrans(gradP[i], wk1[i]);
                    }
                    else
                    {
                        Vmath::Vcopy(physTot, gradP[i], 1, wk1[i], 1);
                    }
                }
                m_mapping->RaiseIndex(wk1, wk2);   // G(p)
           
                m_mapping->Divergence(wk2, F_corrected);   // div(G(p))
                if (!m_mapping->HasConstantJacobian())
                {
Douglas Serson's avatar
Douglas Serson committed
510 511 512
                    Vmath::Vmul(physTot, F_corrected, 1, 
                                         Jac, 1, 
                                         F_corrected, 1);
513
                }                
Douglas Serson's avatar
Douglas Serson committed
514
                // alpha*J*div(G(p))
Douglas Serson's avatar
Douglas Serson committed
515 516
                Vmath::Smul(physTot, m_pressureRelaxation, F_corrected, 1, 
                                                           F_corrected, 1); 
517 518 519 520 521 522 523
                if(m_pressure->GetWaveSpace())
                {
                    m_pressure->HomogeneousFwdTrans(F_corrected, F_corrected);
                }            
                // alpha*J*div(G(p)) - p_ii      
                for (int i = 0; i < m_nConvectiveFields; ++i)
                {
Douglas Serson's avatar
Douglas Serson committed
524 525 526 527
                    m_pressure->PhysDeriv(MultiRegions::DirCartesianMap[i],
                                                            gradP[i], wk1[0]);
                    Vmath::Vsub(physTot, F_corrected, 1, wk1[0], 1, 
                                                            F_corrected, 1);
528 529 530 531
                }
                // p_i,i - J*div(G(p))
                Vmath::Neg(physTot, F_corrected, 1);           
                // alpha*F -  alpha*J*div(G(p)) + p_i,i
Douglas Serson's avatar
Douglas Serson committed
532 533
                Vmath::Smul(physTot, m_pressureRelaxation, Forcing, 1, 
                                                            wk1[0], 1);
534 535 536 537 538
                Vmath::Vadd(physTot, wk1[0], 1, F_corrected, 1, F_corrected, 1);

                //
                // Solve system
                //
Douglas Serson's avatar
Douglas Serson committed
539 540 541 542
                m_pressure->HelmSolve(F_corrected, m_pressure->UpdateCoeffs(), 
                                        NullFlagList, factors);  
                m_pressure->BwdTrans(m_pressure->GetCoeffs(),
                                    m_pressure->UpdatePhys());     
543

544 545 546 547 548 549
                //
                // Test convergence
                //
                error = m_pressure->L2(m_pressure->GetPhys(), previous_iter);           
                if ( forcing_L2 != 0)
                {
550
                    if ( (error/forcing_L2 < m_pressureTolerance))
551 552 553 554 555 556
                    {
                        converged = true;
                    }  
                }
                else
                {
557
                    if ( error < m_pressureTolerance)
558 559 560 561 562
                    {
                        converged = true;
                    }                
                }                    
            }
Douglas Serson's avatar
Douglas Serson committed
563
            if (m_verbose && m_session->GetComm()->GetRank()==0)
564
            {
Douglas Serson's avatar
Douglas Serson committed
565 566
                std::cout << " Pressure system (mapping) converged in " << s <<
                            " iterations with error = " << error << std::endl;
567 568
            }            
        }
569 570 571 572 573 574 575 576 577 578
    }
    
    /**
     * Solve velocity system
     */
    void   VCSMapping::v_SolveViscous(
        const Array<OneD, const Array<OneD, NekDouble> > &Forcing,
        Array<OneD, Array<OneD, NekDouble> > &outarray,
        const NekDouble aii_Dt)
    {
579
        if(!m_implicitViscous)
580
        {
581
            VelocityCorrectionScheme::v_SolveViscous(Forcing, outarray, aii_Dt);
582
        }
583
        else
584
        {
585 586
            int physTot = m_fields[0]->GetTotPoints();
            int nvel = m_nConvectiveFields;
Douglas Serson's avatar
Douglas Serson committed
587 588 589
            bool converged = false;          // flag to mark if system converged
            int s = 0;                       // iteration counter
            NekDouble error, max_error;      // L2 error at current iteration
590
            
Douglas Serson's avatar
Douglas Serson committed
591 592 593
            int maxIter;
            m_session->LoadParameter("MappingMaxIter",maxIter,5000);
            
Douglas Serson's avatar
Douglas Serson committed
594 595
            //L2 norm of F
            Array<OneD, NekDouble> forcing_L2(m_nConvectiveFields,0.0); 
596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611

            // rhs of the equation at current iteration
            Array<OneD, Array<OneD, NekDouble> > F_corrected(nvel);
            // Solution at previous iteration
            Array<OneD, Array<OneD, NekDouble> > previous_iter(nvel);
            // Working space
            Array<OneD, Array<OneD, NekDouble> >  wk(nvel);
            for(int i = 0; i < nvel; ++i)
            {
                F_corrected[i]   = Array<OneD, NekDouble> (physTot, 0.0);
                previous_iter[i] = Array<OneD, NekDouble> (physTot, 0.0);
                wk[i]            = Array<OneD, NekDouble> (physTot, 0.0);
            }

            // Factors for Helmholtz system
            StdRegions::ConstFactorMap factors;
Douglas Serson's avatar
Douglas Serson committed
612 613
            factors[StdRegions::eFactorLambda] = 
                                    1.0*m_viscousRelaxation/aii_Dt/m_kinvis;
614 615 616
            if(m_useSpecVanVisc)
            {
                factors[StdRegions::eFactorSVVCutoffRatio] = m_sVVCutoffRatio;
Douglas Serson's avatar
Douglas Serson committed
617 618
                factors[StdRegions::eFactorSVVDiffCoeff]   = 
                                                      m_sVVDiffCoeff/m_kinvis;
619 620 621 622 623 624
            }

            // Calculate L2-norm of F and set initial solution for iteration
            for(int i = 0; i < nvel; ++i)
            {
                forcing_L2[i] = m_fields[0]->L2(Forcing[i],wk[0]);
Douglas Serson's avatar
Douglas Serson committed
625 626
                m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
                                        previous_iter[i]);
627 628 629 630 631 632 633
            }

            while (!converged)
            {
                converged = true;
                // Iteration counter
                s++;
Douglas Serson's avatar
Douglas Serson committed
634 635 636
                ASSERTL0(s < maxIter,
                         "VCSMapping exceeded maximum number of iterations.");
                
637 638 639 640 641 642
                max_error = 0.0;

                //
                // Calculate forcing term for next iteration
                //

Douglas Serson's avatar
Douglas Serson committed
643
                // Calculate L(U)- in this parts all components might be coupled
644 645 646 647
                if(m_fields[0]->GetWaveSpace())
                {
                    for (int i = 0; i < nvel; ++i)
                    {
Douglas Serson's avatar
Douglas Serson committed
648 649
                        m_fields[0]->HomogeneousBwdTrans(previous_iter[i], 
                                                                    wk[i]);
650 651 652 653 654 655 656 657 658 659
                    }                    
                }
                else
                {
                    for (int i = 0; i < nvel; ++i)
                    {
                        Vmath::Vcopy(physTot, previous_iter[i], 1, wk[i], 1);
                    }                   
                }                    
                
660 661 662
                // (L(U^i) - 1/alpha*U^i_jj)
                m_mapping->VelocityLaplacian(wk, F_corrected,
                                                1.0/m_viscousRelaxation);
663 664 665 666 667
                
                if(m_fields[0]->GetWaveSpace())
                {
                    for (int i = 0; i < nvel; ++i)
                    {
Douglas Serson's avatar
Douglas Serson committed
668 669
                        m_fields[0]->HomogeneousFwdTrans(F_corrected[i], 
                                                            F_corrected[i]);
670 671 672 673 674 675
                    }                    
                }
                else
                {
                    for (int i = 0; i < nvel; ++i)
                    {
Douglas Serson's avatar
Douglas Serson committed
676 677
                        Vmath::Vcopy(physTot, F_corrected[i], 1, 
                                                F_corrected[i], 1);
678 679 680 681 682 683 684
                    }                   
                }
                
                // Loop velocity components
                for (int i = 0; i < nvel; ++i)
                {
                    // (-alpha*L(U^i) + U^i_jj)
685 686 687
                    Vmath::Smul(physTot, -1.0*m_viscousRelaxation,
                                                    F_corrected[i], 1,
                                                    F_corrected[i], 1);
688
                    //  F_corrected = alpha*F + (-alpha*L(U^i) + U^i_jj)
Douglas Serson's avatar
Douglas Serson committed
689 690
                    Vmath::Smul(physTot, m_viscousRelaxation, Forcing[i], 1, 
                                                                    wk[0], 1);
Douglas Serson's avatar
Douglas Serson committed
691
                    Vmath::Vadd(physTot, wk[0], 1, F_corrected[i], 1, 
Douglas Serson's avatar
Douglas Serson committed
692
                                                    F_corrected[i], 1);
693 694 695 696

                    //
                    // Solve System
                    //
Douglas Serson's avatar
Douglas Serson committed
697 698 699
                    m_fields[i]->HelmSolve(F_corrected[i], 
                                    m_fields[i]->UpdateCoeffs(),
                                    NullFlagList, factors); 
700 701 702 703 704 705 706 707 708
                    m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),outarray[i]);

                    //
                    // Test convergence
                    //
                    error = m_fields[i]->L2(outarray[i], previous_iter[i]);

                    if ( forcing_L2[i] != 0)
                    {
709
                        if ( (error/forcing_L2[i] >= m_viscousTolerance))
710 711
                        {
                            converged = false;
Douglas Serson's avatar
Douglas Serson committed
712
                        }
713 714 715
                    }
                    else
                    {
716
                        if ( error >= m_viscousTolerance)
717 718
                        {
                            converged = false;
Douglas Serson's avatar
Douglas Serson committed
719 720
                        }
                    }
721 722 723 724 725 726 727 728 729
                    if (error > max_error)
                    {
                        max_error = error;
                    }

                    // Copy field to previous_iter
                    Vmath::Vcopy(physTot, outarray[i], 1, previous_iter[i], 1); 
                }           
            }
Douglas Serson's avatar
Douglas Serson committed
730
            if (m_verbose && m_session->GetComm()->GetRank()==0)
731
            {
Douglas Serson's avatar
Douglas Serson committed
732 733
                std::cout << " Velocity system (mapping) converged in " << s <<
                        " iterations with error = " << max_error << std::endl;
734
            }            
735 736 737 738 739 740 741
        }
    }
    
    /**
     * Explicit terms of the mapping
     */
    void   VCSMapping::ApplyIncNSMappingForcing(
742
        const Array<OneD, Array<OneD, NekDouble> >    &inarray,
743 744 745
        Array<OneD, Array<OneD, NekDouble> >          &outarray)
    {
        int physTot = m_fields[0]->GetTotPoints();
Douglas Serson's avatar
Douglas Serson committed
746 747 748 749
        Array<OneD, Array<OneD, NekDouble> >       vel(m_nConvectiveFields);
        Array<OneD, Array<OneD, NekDouble> >       velPhys(m_nConvectiveFields);
        Array<OneD, Array<OneD, NekDouble> >       Forcing(m_nConvectiveFields);
        Array<OneD, Array<OneD, NekDouble> >       tmp(m_nConvectiveFields);
750 751
        for (int i = 0; i < m_nConvectiveFields; ++i)
        {
752
            velPhys[i] = Array<OneD, NekDouble> (physTot, 0.0);
753 754 755 756
            Forcing[i] = Array<OneD, NekDouble> (physTot, 0.0);
            tmp[i] = Array<OneD, NekDouble> (physTot, 0.0);
        }
        
757
        // Get fields and store velocity in wavespace and physical space
758 759 760 761
        if(m_fields[0]->GetWaveSpace())
        {
            for (int i = 0; i < m_nConvectiveFields; ++i)
            {
762
                vel[i] = inarray[i];
763
                m_fields[0]->HomogeneousBwdTrans(vel[i],velPhys[i]);
764 765 766 767 768 769
            }
        }
        else
        {
            for (int i = 0; i < m_nConvectiveFields; ++i)
            {
770 771
                vel[i] = inarray[i];
                Vmath::Vcopy(physTot, inarray[i], 1, velPhys[i], 1);
772 773
            }
        }
774

775
        //Advection contribution
776
        MappingAdvectionCorrection(velPhys, Forcing);
777 778 779 780
        
        // Time-derivative contribution
        if ( m_mapping->IsTimeDependent() )
        {
781
            MappingAccelerationCorrection(vel, velPhys, tmp);
782 783 784 785 786 787
            for (int i = 0; i < m_nConvectiveFields; ++i)
            {
                Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
            }            
        }
        
788
        // Pressure contribution
789
        if (!m_implicitPressure)
790
        {
791
            MappingPressureCorrection(tmp);
792 793 794 795
            for (int i = 0; i < m_nConvectiveFields; ++i)
            {
                Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
            }            
796 797
        }   
        // Viscous contribution
798
        if ( (!m_implicitViscous) && (!m_neglectViscous))
799
        {
800
            MappingViscousCorrection(velPhys, tmp);
801 802 803 804 805 806
            for (int i = 0; i < m_nConvectiveFields; ++i)
            {            
                Vmath::Smul(physTot, m_kinvis, tmp[i], 1, tmp[i], 1);
                Vmath::Vadd(physTot, tmp[i], 1, Forcing[i], 1, Forcing[i], 1);
            }
        }        
807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823

        // If necessary, transform to wavespace
        if(m_fields[0]->GetWaveSpace())
        {
            for (int i = 0; i < m_nConvectiveFields; ++i)
            {
                m_fields[0]->HomogeneousFwdTrans(Forcing[i],Forcing[i]);
            }
        }
        
        // Add to outarray 
        for (int i = 0; i < m_nConvectiveFields; ++i)
        {
            Vmath::Vadd(physTot, outarray[i], 1, Forcing[i], 1, outarray[i], 1);
        }          
    }
    
824
        void VCSMapping::MappingAdvectionCorrection(
825
            const Array<OneD, Array<OneD, NekDouble> >        &velPhys,
826 827 828 829 830 831 832 833
            Array<OneD, Array<OneD, NekDouble> >              &outarray)
        {
            int physTot = m_fields[0]->GetTotPoints();
            int nvel = m_nConvectiveFields;
            
            Array<OneD, Array<OneD, NekDouble> > wk(nvel*nvel);

            // Apply Christoffel symbols to obtain {i,kj}vel(k) 
834
            m_mapping->ApplyChristoffelContravar(velPhys, wk);
835 836 837 838 839 840 841

            // Calculate correction -U^j*{i,kj}vel(k)
            for (int i = 0; i< nvel; i++)
            {
                Vmath::Zero(physTot,outarray[i],1);
                for (int j = 0; j< nvel; j++)
                {
842
                        Vmath::Vvtvp(physTot,wk[i*nvel+j],1,velPhys[j],1,
843 844 845 846 847 848 849
                                        outarray[i],1,outarray[i],1);
                }    
                Vmath::Neg(physTot, outarray[i], 1);
            } 
        }
        
        void VCSMapping::MappingAccelerationCorrection(
850 851
            const Array<OneD, Array<OneD, NekDouble> >        &vel,
            const Array<OneD, Array<OneD, NekDouble> >        &velPhys,
852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869
            Array<OneD, Array<OneD, NekDouble> >              &outarray)
        {
            int physTot = m_fields[0]->GetTotPoints();
            int nvel = m_nConvectiveFields;
            
            Array<OneD, Array<OneD, NekDouble> > wk(nvel*nvel);
            Array<OneD, Array<OneD, NekDouble> > tmp(nvel);
            Array<OneD, Array<OneD, NekDouble> > coordVel(nvel);
            for (int i = 0; i< nvel; i++)
            {
                tmp[i] = Array<OneD, NekDouble> (physTot, 0.0);
                coordVel[i] = Array<OneD, NekDouble> (physTot, 0.0);
            }
            // Get coordinates velocity in transformed system
            m_mapping->GetCoordVelocity(tmp);
            m_mapping->ContravarFromCartesian(tmp, coordVel);         
            
            // Calculate first term: U^j u^i,j = U^j (du^i/dx^j + {i,kj}u^k)
870
            m_mapping->ApplyChristoffelContravar(velPhys, wk);
871 872 873
            for (int i=0; i< nvel; i++)
            {
                Vmath::Zero(physTot,outarray[i],1);
874 875

                m_fields[0]->PhysDeriv(velPhys[i], tmp[0], tmp[1]);
876 877
                for (int j=0; j< nvel; j++)
                {
878
                    if (j == 2)
879 880
                    {
                        m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j],
881
                                        vel[i], tmp[2]);
882 883
                        if (m_fields[0]->GetWaveSpace())
                        {
884
                            m_fields[0]->HomogeneousBwdTrans(tmp[2],tmp[2]);
885 886
                        }
                    }
887

888
                    Vmath::Vadd(physTot,wk[i*nvel+j],1,tmp[j],1,
Douglas Serson's avatar
Douglas Serson committed
889
                                                        wk[i*nvel+j], 1); 
890 891 892 893 894 895
                    
                    Vmath::Vvtvp(physTot, coordVel[j], 1, wk[i*nvel+j], 1,
                                          outarray[i], 1, outarray[i], 1);
                }
            }
            
896 897 898 899
            // Set wavespace to false and store current value
            bool wavespace = m_fields[0]->GetWaveSpace();
            m_fields[0]->SetWaveSpace(false);

900
            // Add -u^j U^i,j
901
            m_mapping->ApplyChristoffelContravar(coordVel, wk);
902 903
            for (int i=0; i< nvel; i++)
            {
904
                if(nvel == 2)
905
                {
906 907 908 909 910 911
                    m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1]);
                }
                else
                {
                    m_fields[0]->PhysDeriv(coordVel[i], tmp[0], tmp[1], tmp[2]);
                }
912

913 914 915
                for (int j=0; j< nvel; j++)
                {
                    Vmath::Vadd(physTot,wk[i*nvel+j],1,tmp[j],1,
Douglas Serson's avatar
Douglas Serson committed
916
                                                        wk[i*nvel+j], 1);
917 918
                    Vmath::Neg(physTot, wk[i*nvel+j], 1);
                    
919
                    Vmath::Vvtvp(physTot, velPhys[j], 1, wk[i*nvel+j], 1,
920 921 922
                                          outarray[i], 1, outarray[i], 1);
                }
            }
Douglas Serson's avatar
Douglas Serson committed
923

924 925 926 927 928 929 930 931 932 933
            // Restore value of wavespace 
            m_fields[0]->SetWaveSpace(wavespace);
        }

        void VCSMapping::MappingPressureCorrection(
            Array<OneD, Array<OneD, NekDouble> >              &outarray)
        {
            int physTot = m_fields[0]->GetTotPoints();
            int nvel = m_nConvectiveFields;

934 935
            // Calculate g^(ij)p_(,j)
            m_mapping->RaiseIndex(m_gradP, outarray);
Douglas Serson's avatar
Douglas Serson committed
936

937 938 939 940 941 942 943 944
            // Calculate correction = (nabla p)/J - g^(ij)p_,j
            // (Jac is not required if it is constant)
            if ( !m_mapping->HasConstantJacobian())
            {
                Array<OneD, NekDouble> Jac(physTot, 0.0);
                m_mapping->GetJacobian(Jac);
                for(int i = 0; i < nvel; ++i)
                {
945
                    Vmath::Vdiv(physTot, m_gradP[i], 1, Jac, 1, m_gradP[i], 1);
946 947 948 949
                }
            }
            for(int i = 0; i < nvel; ++i)
            {
950 951
                Vmath::Vsub(physTot, m_gradP[i], 1,outarray[i], 1,
                                                    outarray[i],1);
952 953 954 955
            }
        }

        void VCSMapping::MappingViscousCorrection(
956
            const Array<OneD, Array<OneD, NekDouble> >        &velPhys,
957 958
            Array<OneD, Array<OneD, NekDouble> >              &outarray)
        {
959 960
            // L(U) - 1.0*d^2(u^i)/dx^jdx^j
            m_mapping->VelocityLaplacian(velPhys, outarray, 1.0);
961
        }
Douglas Serson's avatar
Douglas Serson committed
962

963
} //end of namespace