Commit aef14d80 authored by Spencer Sherwin's avatar Spencer Sherwin

Added routine to calculate pressure normalisation based on Linf rather than L2 norm


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@3775 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 5338b173
......@@ -18,6 +18,26 @@ SET(VortexWaveInteractionSolverSource
./VortexWaveInteractionSolver.cpp
)
SET(CalcL2ToLinfPressureSource
../Auxiliary/Driver.cpp
../Auxiliary/DriverStandard.cpp
../Auxiliary/DriverArnoldi.cpp
../Auxiliary/DriverModifiedArnoldi.cpp
../Auxiliary/EquationSystem.cpp
../ADRSolver/EquationSystems/SteadyAdvectionDiffusion.cpp
../ADRSolver/EquationSystems/SteadyAdvectionDiffusionReaction.cpp
../IncNavierStokesSolver/EquationSystems/CoupledLinearNS.cpp
../IncNavierStokesSolver/EquationSystems/CoupledLocalToGlobalC0ContMap.cpp
../IncNavierStokesSolver/EquationSystems/IncNavierStokes.cpp
../IncNavierStokesSolver/EquationSystems/VelocityCorrectionScheme.cpp
../IncNavierStokesSolver/AdvectionTerms/AdvectionTerm.cpp
../IncNavierStokesSolver/AdvectionTerms/NavierStokesAdvection.cpp
../IncNavierStokesSolver/AdvectionTerms/LinearisedAdvection.cpp
../IncNavierStokesSolver/AdvectionTerms/AdjointAdvection.cpp
./CalcL2toLinfPressure.cpp
./VortexWaveInteraction.cpp
)
IF (NEKTAR_USE_ARPACK)
SET(VortexWaveInteractionSolverSource ${VortexWaveInteractionSolverSource}
../Auxiliary/DriverArpack.cpp)
......@@ -25,4 +45,7 @@ ENDIF (NEKTAR_USE_ARPACK)
ADD_SOLVER_EXECUTABLE(VortexWaveInteractionSolver solvers
${VortexWaveInteractionSolverSource})
ADD_SOLVER_EXECUTABLE(CalcL2ToLinfPressure solvers
${CalcL2ToLinfPressureSource})
///////////////////////////////////////////////////////////////////////////////
//
// File CalcL2ToLinfPressure.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: Vortex Wave Interaction conversion from L2 to Linf
// Pressure Normalisation
//
///////////////////////////////////////////////////////////////////////////////
#include <Auxiliary/Driver.h>
#include <LibUtilities/BasicUtils/SessionReader.h>
using namespace Nektar;
#include<VortexWaveInteraction/VortexWaveInteraction.h>
int main(int argc, char *argv[])
{
try
{
VortexWaveInteraction vwi(argc,argv);
vwi.CalcL2ToLinfPressure();
return 1;
}
catch (const std::runtime_error&)
{
return 1;
}
catch (const std::string& eStr)
{
cout << "Error: " << eStr << endl;
}
}
......@@ -83,6 +83,15 @@ namespace Nektar
m_deltaFcnDecay = 0.0;
}
if(m_sessionVWI->DefinesSolverInfo("LinfPressureNorm"))
{
m_useLinfPressureNorm = true;
}
else
{
m_useLinfPressureNorm = false;
}
if( m_sessionVWI->DefinesSolverInfo("INTERFACE") )
{
m_iterinterface = true;
......@@ -93,8 +102,8 @@ namespace Nektar
}
m_alpha = Array<OneD, NekDouble> (storesize);
m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
m_waveForceMag = Array<OneD, NekDouble> (storesize);
m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
m_waveForceMag = Array<OneD, NekDouble> (storesize);
m_waveForceMag[0] = m_sessionVWI->GetParameter("WaveForceMag");
m_leading_real_evl = Array<OneD, NekDouble> (storesize);
......@@ -554,29 +563,41 @@ namespace Nektar
m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),
m_wavePressure->GetPlane(1)->UpdatePhys());
// Determine normalisation of pressure so that |P|/A = 1
NekDouble norm = 0, l2;
l2 = m_wavePressure->GetPlane(0)->L2();
norm = l2*l2;
l2 = m_wavePressure->GetPlane(1)->L2();
norm += l2*l2;
Vmath::Fill(2*npts,1.0,der1,1);
NekDouble area = m_waveVelocities[0]->GetPlane(0)->PhysIntegral(der1);
norm = sqrt(area/norm);
NekDouble invnorm;
if(m_useLinfPressureNorm)
{
NekDouble Linf;
Linf = m_wavePressure->Linf(m_wavePressure->GetPhys());
invnorm = 1.0/Linf;
}
else
{
// Determine normalisation of pressure so that |P|/A = 1
NekDouble l2;
l2 = m_wavePressure->GetPlane(0)->L2();
invnorm = l2*l2;
l2 = m_wavePressure->GetPlane(1)->L2();
invnorm += l2*l2;
Vmath::Fill(2*npts,1.0,der1,1);
NekDouble area = m_waveVelocities[0]->GetPlane(0)->PhysIntegral(der1);
invnorm = sqrt(area/invnorm);
}
// Get hold of arrays.
m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
Array<OneD, NekDouble> u_real = m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
Vmath::Smul(npts,norm,u_real,1,u_real,1);
Vmath::Smul(npts,invnorm,u_real,1,u_real,1);
m_waveVelocities[0]->GetPlane(1)->BwdTrans(m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
Array<OneD, NekDouble> u_imag = m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
Vmath::Smul(npts,norm,u_imag,1,u_imag,1);
Vmath::Smul(npts,invnorm,u_imag,1,u_imag,1);
m_waveVelocities[1]->GetPlane(0)->BwdTrans(m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
Array<OneD, NekDouble> v_real = m_waveVelocities[1]->GetPlane(0)->UpdatePhys();
Vmath::Smul(npts,norm,v_real,1,v_real,1);
Vmath::Smul(npts,invnorm,v_real,1,v_real,1);
m_waveVelocities[1]->GetPlane(1)->BwdTrans(m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
Array<OneD, NekDouble> v_imag = m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
Vmath::Smul(npts,norm,v_imag,1,v_imag,1);
Vmath::Smul(npts,invnorm,v_imag,1,v_imag,1);
// Calculate non-linear terms for x and y directions
// d/dx(u u* + u* u)
......@@ -742,6 +763,33 @@ namespace Nektar
}
}
void VortexWaveInteraction::CalcL2ToLinfPressure(void)
{
ExecuteWave();
NekDouble Linf;
Linf = m_wavePressure->Linf(m_wavePressure->GetPhys());
cout << "Linf: " << Linf << endl;
NekDouble l2,norm;
l2 = m_wavePressure->GetPlane(0)->L2();
norm = l2*l2;
l2 = m_wavePressure->GetPlane(1)->L2();
norm += l2*l2;
int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
Array<OneD, NekDouble> der1(npts);
Vmath::Fill(2*npts,1.0,der1,1);
NekDouble area = m_waveVelocities[0]->GetPlane(0)->PhysIntegral(der1);
l2 = sqrt(norm/area);
cout << "L2: " << l2 << endl;
cout << "Ratio Linf/L2: "<< Linf/l2 << endl;
}
void VortexWaveInteraction::SaveFile(string file, string dir, int n)
{
static map<string,int> opendir;
......@@ -1561,8 +1609,5 @@ cout<<"ucnt="<<cnt<<endl;
cnt++;
}
}
}
......@@ -77,6 +77,8 @@ namespace Nektar
void SaveLoopDetails(string dir, int i);
void CalcNonLinearWaveForce(void);
void CalcL2ToLinfPressure (void);
void SaveFile(string fileend, string dir, int n);
void MoveFile(string fileend, string dir, int n);
void CopyFile(string file1end, string file2end);
......@@ -175,6 +177,8 @@ namespace Nektar
int m_maxWaveForceMagIter;
bool m_deltaFcnApprox; // Activate delta function approximation around wave
bool m_useLinfPressureNorm; // Activate if use Pressure Linf Normalisation
NekDouble m_deltaFcnDecay; // Delta function decay level
Array<OneD, NekDouble> m_waveForceMag;
......
......@@ -53,7 +53,7 @@ int main(int argc, char *argv[])
try
{
VortexWaveInteraction vwi(argc,argv);
for(int i = 0; i < vwi.GetMaxWaveForceMagIter(); ++i)
{
DoFixedForcingIteration(vwi);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment