VortexWaveInteractionSolver.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: VortexWaveInteractionSolver.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).
//
// 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: Vortex Wave Interaction solver
//
///////////////////////////////////////////////////////////////////////////////
 
#include <VortexWaveInteraction/VortexWaveInteraction.h>
 
using namespace std;
using namespace Nektar;
 
void DoFixedForcingIteration(VortexWaveInteraction &vwi);
void Mvdir(string dir, NekDouble dir_ending);
 
int main(int argc, char *argv[])
{
    try
    {
        VortexWaveInteraction vwi(argc, argv);
 
        for (int i = 0; i < vwi.GetMaxWaveForceMagIter(); ++i)
        {
            DoFixedForcingIteration(vwi);
 
            NekDouble WaveForce = vwi.GetWaveForceMag();
 
            vwi.AppendEvlToFile("ConvergedSolns", WaveForce);
 
            vwi.SetWaveForceMag(WaveForce + vwi.GetWaveForceMagStep());
            vwi.SetPrevAlpha(vwi.GetAlpha());
            // vwi.SetAlpha(vwi.GetAlpha() +
            // vwi.GetDAlphaDWaveForceMag()*vwi.GetWaveForceMagStep());
            vwi.SetAlpha(vwi.GetAlpha() + (vwi.GetWaveForceMagStep() > 0
                                               ? vwi.GetAlphaStep()
                                               : (-vwi.GetAlphaStep())));
 
            // Save data directories.
            if (vwi.GetVWIIterationType() == eFixedAlphaWaveForcing)
            {
                Mvdir("Save", WaveForce);
            }
            else
            {
                Mvdir("Save_Outer", WaveForce);
                // Execute Another loop so that not same initial conditions as
                // last iteration
                vwi.ExecuteLoop();
            }
        }
    }
 
    catch (const std::runtime_error &)
    {
        return 1;
    }
 
    catch (const std::string &eStr)
    {
        cout << "Error: " << eStr << endl;
    }
}
 
void Mvdir(string dir, NekDouble dir_ending)
{
    // save OuterIter.his if exists
    string saveOuterIter = "mv -f OuterIter.his " + dir;
    if (system(saveOuterIter.c_str()))
    {
        ASSERTL0(false, saveOuterIter.c_str());
    }
 
    // Mv directory
    string newdir  = dir + boost::lexical_cast<std::string>(dir_ending);
    string syscall = "mv -f " + dir + " " + newdir;
 
    if (system(syscall.c_str()))
    {
        ASSERTL0(false, syscall.c_str());
    }
 
    // make new directory
    syscall = "mkdir " + dir;
    ASSERTL0(system(syscall.c_str()) == 0,
             "Failed to make directory '" + dir + "'");
}
 
void DoFixedForcingIteration(VortexWaveInteraction &vwi)
{
 
    // Reset eigenvalue checker in case used in previous iterations
    vwi.CheckEigIsStationary(true);
 
    switch (vwi.GetVWIIterationType())
    {
        case eFixedAlpha:
        {
            int i;
            int nouter_iter     = vwi.GetNOuterIterations();
            bool exit_iteration = false;
 
            while (exit_iteration == false)
            {
                // Reset eigenvalue checker in case used in previous iterations
                vwi.CheckEigIsStationary(true);
 
                for (i = vwi.GetIterStart(); i < vwi.GetIterEnd(); ++i)
                {
                    vwi.ExecuteLoop();
                    vwi.SaveLoopDetails("Save", i);
                    vwi.AppendEvlToFile("conv.his", i);
 
                    if (vwi.CheckEigIsStationary())
                    {
                        vwi.SaveLoopDetails("Save_Outer", nouter_iter);
                        break;
                    }
                }
 
                // check to see if growth was converged.
                if (i == vwi.GetIterEnd())
                {
                    cout << "Failed to converge growth rate in"
                         << " inner iteration after " << vwi.GetIterEnd()
                         << " loops" << endl;
                    exit(1);
                }
 
                vwi.AppendEvlToFile("OuterIter.his", nouter_iter++);
                exit_iteration = vwi.CheckIfAtNeutralPoint();
                if (exit_iteration == false)
                {
                    vwi.UpdateWaveForceMag(nouter_iter);
                }
 
                if (nouter_iter >= vwi.GetMaxOuterIterations())
                {
                    cerr << "Failed to converge after "
                         << vwi.GetMaxOuterIterations() << " outer iterations"
                         << endl;
                    exit_iteration = true;
                }
            }
        }
        break;
        case eFixedWaveForcing:
        {
            int i;
            int nouter_iter         = vwi.GetNOuterIterations();
            bool exit_iteration     = false;
            NekDouble saveEigRelTol = vwi.GetEigRelTol();
            int saveMinIters        = vwi.GetMinInnerIterations();
            int init_search         = 0;
 
            if (init_search)
            {
                // initial set m_eigelTol to 1e-1 and inner iterations to 1 for
                // quick search
                vwi.SetEigRelTol(1e-1);
                vwi.SetMinInnerIterations(2);
            }
 
            while (exit_iteration == false)
            {
                // Reset eigenvalue checker in case used in previous iterations
                vwi.CheckEigIsStationary(true);
 
                for (i = vwi.GetIterStart(); i < vwi.GetIterEnd(); ++i)
                {
                    vwi.ExecuteLoop();
                    vwi.SaveLoopDetails("Save", i);
                    vwi.AppendEvlToFile("conv.his", i);
 
                    if (vwi.CheckEigIsStationary())
                    {
                        vwi.SaveLoopDetails("Save_Outer", nouter_iter);
                        break;
                    }
                }
 
                // check to see if growth was converged.
                if (i == vwi.GetIterEnd())
                {
                    cout << "Failed to converge growth rate in"
                         << " inner iteration after " << vwi.GetIterEnd()
                         << " loops" << endl;
                    exit(1);
                }
 
                vwi.AppendEvlToFile("OuterIter.his", nouter_iter++);
 
                exit_iteration = vwi.CheckIfAtNeutralPoint();
 
                cout << "m_alpha[0] is " << vwi.GetAlpha() << endl;
 
                if (exit_iteration == false)
                {
                    vwi.UpdateAlpha(nouter_iter);
                    if (vwi.IfIterInterface() == true)
                    {
                        vwi.CalcNonLinearWaveForce();
                    }
                }
 
                // Redo iteration if at first coarse search
                if ((exit_iteration == true) && (init_search == 1))
                {
                    init_search = 0;
                    vwi.SetEigRelTol(saveEigRelTol);
                    vwi.SetMinInnerIterations(saveMinIters);
                    nouter_iter    = 1;
                    exit_iteration = false;
                }
 
                if (nouter_iter >= vwi.GetMaxOuterIterations())
                {
                    cerr << "Failed to converge after "
                         << vwi.GetMaxOuterIterations() << " outer iterations"
                         << endl;
                    exit_iteration = true;
                }
            }
 
            vwi.UpdateDAlphaDWaveForceMag(vwi.GetPrevAlpha());
        }
        break;
        case eFixedAlphaWaveForcing:
 
            for (int i = vwi.GetIterStart(); i < vwi.GetIterEnd(); ++i)
            {
                vwi.ExecuteLoop();
                vwi.SaveLoopDetails("Save", i);
                vwi.AppendEvlToFile("conv.his", i);
            }
            break;
        case eFixedWaveForcingWithSubIterationOnAlpha:
        {
            int i;
            int nouter_iter     = vwi.GetNOuterIterations();
            bool exit_iteration = false;
 
            while (exit_iteration == false)
            {
                bool exit_alphaIter = false;
 
                vwi.ExecuteLoop(false);
 
                // Sub iterate Alpha
                for (i = 0; i < vwi.GetIterEnd(); ++i)
                {
                    vwi.SaveLoopDetails("Save", i);
 
                    vwi.AppendEvlToFile("AlphaIter.his", i);
 
                    exit_alphaIter = vwi.CheckIfAtNeutralPoint();
                    if (exit_alphaIter == false)
                    {
                        vwi.UpdateAlpha(i + 1);
                        vwi.ExecuteWave();
                    }
                    else
                    {
                        vwi.CalcNonLinearWaveForce();
                        break;
                    }
                }
 
                // check to see if growth was converged.
                if (i == vwi.GetIterEnd())
                {
                    cout << "Failed to converge growth rate in"
                         << " inner iteration after " << vwi.GetIterEnd()
                         << " loops" << endl;
                    exit(1);
                }
 
                vwi.MoveFile("AlphaIter.his", "Save_Outer", nouter_iter);
                vwi.SaveLoopDetails("Save_Outer", nouter_iter);
                vwi.AppendEvlToFile("OuterIter.his", nouter_iter++);
 
                // assume that if only previous inner loop has
                // only done one iteration then we are at neutral
                // point
                if (i == 0 && vwi.IfIterInterface() == false)
                {
                    exit_iteration = true;
                }
 
                if (nouter_iter >= vwi.GetMaxOuterIterations())
                {
                    cerr << "Failed to converge after "
                         << vwi.GetMaxOuterIterations() << " outer iterations"
                         << endl;
                    exit_iteration = true;
                }
            }
        }
        break;
        default:
            ASSERTL0(false, "Unknown iteration type");
    }
}