Commit 82bcac00 authored by Robin Basso's avatar Robin Basso

test arnoldi adjoint

parent 3515af31
......@@ -84,10 +84,10 @@ void DriverArnoldi::v_InitObject(ostream &out)
m_BlowingSuction = m_equ[0]->v_CheckBSBC();
// Check to see if a BS-BC is defined
if(m_BlowingSuction == true)
if(m_BlowingSuction == true && m_session->GetSolverInfo("EvolutionOperator") == "Direct")
{
// Determine number of BS-BC fields the solve is supporting;
m_nBSBCFields = 2;
m_nBSBCFields = 0;
}
else
{
......@@ -151,50 +151,51 @@ void DriverArnoldi::v_InitObject(ostream &out)
boost::iequals(m_session->GetSolverInfo("ModeType"),
"SingleMode"))
{
out << "\tSingle Fourier mode : true " << endl;
out << "\tSingle Fourier mode : true " << endl;
ASSERTL0(m_session->DefinesSolverInfo("Homogeneous"),
"Expected a homogeneous expansion to be defined "
"with single mode");
}
else
{
out << "\tSingle Fourier mode : false " << endl;
out << "\tSingle Fourier mode : false " << endl;
}
if(m_session->DefinesSolverInfo("BetaZero"))
{
out << "\tBeta set to Zero : true (overrides LHom)"
out << "\tBeta set to Zero : true (overrides LHom)"
<< endl;
}
else
{
out << "\tBeta set to Zero : false " << endl;
out << "\tBeta set to Zero : false " << endl;
}
if(m_timeSteppingAlgorithm)
{
out << "\tEvolution operator : "
out << "\tEvolution operator : "
<< m_session->GetSolverInfo("EvolutionOperator")
<< endl;
}
else
{
out << "\tShift (Real,Imag) : " << m_realShift
out << "\tShift (Real,Imag) : " << m_realShift
<< "," << m_imagShift << endl;
}
if (m_BlowingSuction == true)
{
out << "\tBlowing/Suction BC's : true" << endl;
out << "\tBlowing/Suction BC's : true" << endl;
}
else
{
out << "\tBlowing/Suction BC's : false" << endl;
out << "\tBlowing/Suction BC's : false" << endl;
}
out << "\tKrylov-space dimension : " << m_kdim << endl;
out << "\tNumber of vectors : " << m_nvec << endl;
out << "\tMax iterations : " << m_nits << endl;
out << "\tEigenvalue tolerance : " << m_evtol << endl;
out << "\tArnoldi iteration's period : " << m_period << endl;
out << "\tKrylov-space dimension : " << m_kdim << endl;
out << "\tNumber of vectors : " << m_nvec << endl;
out << "\tMax iterations : " << m_nits << endl;
out << "\tEigenvalue tolerance : " << m_evtol << endl;
out << "======================================================"
<< endl;
}
......@@ -216,13 +217,13 @@ void DriverArnoldi::v_InitObject(ostream &out)
fields[k]->SetPhysState(false);
}
if(m_BlowingSuction == true)
{
NekDouble theta = array[m_nFields*nq];
NekDouble thetadot = array[(m_nFields+1)*nq];
// if(m_BlowingSuction == true && m_session->GetSolverInfo("EvolutionOperator") == "Direct")
// {
// NekDouble theta = array[m_nFields*nq];
// NekDouble thetadot = array[(m_nFields+1)*nq];
m_equ[0]->v_SetStruct(theta, thetadot);
}
// m_equ[0]->v_SetStruct(theta, thetadot);
// }
};
/**
......@@ -253,17 +254,17 @@ void DriverArnoldi::v_InitObject(ostream &out)
fields[k]->SetPhysState(false);
}
if(m_BlowingSuction == true)
{
NekDouble theta = 0.0;
NekDouble thetadot = 0.0;
m_equ[0]->v_GetStruct(theta, thetadot);
Array<OneD, NekDouble> tmp1(nq, theta);
Array<OneD, NekDouble> tmp2(nq, thetadot);
// if(m_BlowingSuction == true && m_session->GetSolverInfo("EvolutionOperator") == "Direct")
// {
// NekDouble theta = 0.0;
// NekDouble thetadot = 0.0;
// m_equ[0]->v_GetStruct(theta, thetadot);
// Array<OneD, NekDouble> tmp1(nq, theta);
// Array<OneD, NekDouble> tmp2(nq, thetadot);
Vmath::Vcopy(nq, &tmp1[0], 1, &array[m_nFields*nq], 1);
Vmath::Vcopy(nq, &tmp2[0], 1, &array[(m_nFields+1)*nq], 1);
}
// Vmath::Vcopy(nq, &tmp1[0], 1, &array[m_nFields*nq], 1);
// Vmath::Vcopy(nq, &tmp2[0], 1, &array[(m_nFields+1)*nq], 1);
// }
};
......
......@@ -124,6 +124,8 @@ void DriverModifiedArnoldi::v_Execute(ostream &out)
Array<OneD, NekDouble> wr = Array<OneD, NekDouble> (m_kdim, 0.0);
Array<OneD, NekDouble> wi = Array<OneD, NekDouble> (m_kdim, 0.0);
Array<OneD, NekDouble> zvec = Array<OneD, NekDouble> (m_kdim*m_kdim, 0.0);
NekDouble theta = 0.0;
NekDouble thetadot = 0.0;
Array<OneD, Array<OneD, NekDouble> > Kseq
= Array<OneD, Array<OneD, NekDouble> > (m_kdim + 1);
......@@ -166,18 +168,20 @@ void DriverModifiedArnoldi::v_Execute(ostream &out)
}
// Normalise first vector in sequence based on vel field
//alpha[0] = Blas::Ddot(m_nFields*nq, &Kseq[0][0], 1, &Kseq[0][0], 1);
alpha[0] = Blas::Ddot(ntot, &Kseq[0][0], 1, &Kseq[0][0], 1);
m_comm->AllReduce(alpha[0], Nektar::LibUtilities::ReduceSum);
alpha[0] = std::sqrt(alpha[0]);
Vmath::Smul(ntot, 1.0/alpha[0], Kseq[0], 1, Kseq[0], 1);
// if(m_nBSBCFields)
// {
// alpha[0] = Blas::Ddot(ntot, &Kseq[0][0], 1, &Kseq[0][0], 1);
// m_comm->AllReduce(alpha[0], Nektar::LibUtilities::ReduceSum);
// alpha[0] = std::sqrt(alpha[0]);
// }
// Scale angle if needed
if(m_BlowingSuction)
{
m_equ[0]->v_GetStruct(theta, thetadot);
theta = theta/alpha[0];
thetadot = thetadot/alpha[0];
// Send back the angle values to time marching solver
m_equ[0]->v_SetStruct(theta, thetadot);
}
// Fill initial krylov sequence
NekDouble resid0;
......@@ -187,18 +191,20 @@ void DriverModifiedArnoldi::v_Execute(ostream &out)
EV_update(Kseq[i-1], Kseq[i]);
// Normalise based on nfields
//alpha[i] = Blas::Ddot(m_nFields*nq, &Kseq[i][0], 1, &Kseq[i][0], 1);
alpha[i] = Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[i][0], 1);
m_comm->AllReduce(alpha[i], Nektar::LibUtilities::ReduceSum);
alpha[i] = std::sqrt(alpha[i]);
Vmath::Smul(ntot, 1.0/alpha[i], Kseq[i], 1, Kseq[i], 1);
// if(m_nBSBCFields)
// {
// alpha[i] = Blas::Ddot(ntot, &Kseq[i][0], 1, &Kseq[i][0], 1);
// m_comm->AllReduce(alpha[i], Nektar::LibUtilities::ReduceSum);
// alpha[i] = std::sqrt(alpha[i]);
// }
// Scale angle if needed
if(m_BlowingSuction)
{
m_equ[0]->v_GetStruct(theta, thetadot);
theta = theta/alpha[i];
thetadot = thetadot/alpha[i];
// Send back the angle values to time marching solver
m_equ[0]->v_SetStruct(theta, thetadot);
}
// Copy Krylov sequence into temporary storage
for (int k = 0; k < i + 1; ++k)
......@@ -251,6 +257,15 @@ void DriverModifiedArnoldi::v_Execute(ostream &out)
alpha[m_kdim] = std::sqrt(alpha[m_kdim]);
Vmath::Smul(ntot, 1.0/alpha[m_kdim], Kseq[m_kdim], 1,
Kseq[m_kdim], 1);
// Scale angle if needed
if(m_BlowingSuction)
{
m_equ[0]->v_GetStruct(theta, thetadot);
theta = theta/alpha[m_kdim];
thetadot = thetadot/alpha[m_kdim];
// Send back the angle values to time marching solver
m_equ[0]->v_SetStruct(theta, thetadot);
}
// Copy Krylov sequence into temporary storage
for (int k = 0; k < m_kdim + 1; ++k)
......
......@@ -892,12 +892,14 @@ namespace Nektar
////////// TEST angleVel = sigma * exp(sigma * t) ///////
NekDouble sigma = 6.7256e-02;
NekDouble timeCopy = time;
while (timeCopy > m_finTime)
{
timeCopy -= m_finTime;
}
m_bsbcParams->m_angleVel[0] = sigma * exp(sigma * timeCopy);
m_bsbcParams->m_angle = exp(sigma * timeCopy);
m_bsbcParams->m_angleVel[0] = 0.04 * sigma * exp(sigma * timeCopy);
m_bsbcParams->m_angle = exp(sigma * timeCopy);
////////////////////////////////////////////////////////
if( m_bsbcParams->m_doOutput )
......
......@@ -794,13 +794,13 @@ namespace Nektar
// Update BC for blowing/suction simulations:
if (m_BlowingSuction && m_session->GetSolverInfo("EvolutionOperator") == "Direct")
{
ScaleBSBCDirect();
SolveStructuralDirect(m_time);
ScaleBSBCDirect();
}
if (m_BlowingSuction && m_session->GetSolverInfo("EvolutionOperator") == "Adjoint")
{
ScaleBSBCAdjoint();
SolveStructuralAdjoint(m_time);
ScaleBSBCAdjoint();
}
}
......
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