Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Attributes | List of all members
Nektar::VortexWaveInteraction Class Reference

#include <VortexWaveInteraction.h>

Collaboration diagram for Nektar::VortexWaveInteraction:
Collaboration graph
[legend]

Public Member Functions

 VortexWaveInteraction (int argc, char *argv[])
 
 ~VortexWaveInteraction (void)
 
void ExecuteRoll (void)
 
void ExecuteStreak (void)
 
void ExecuteWave (void)
 
void ExecuteLoop (bool CalcWaveForce=true)
 
void SaveLoopDetails (std::string dir, int i)
 
void CalcNonLinearWaveForce (void)
 
void CalcL2ToLinfPressure (void)
 
void SaveFile (std::string fileend, std::string dir, int n)
 
void MoveFile (std::string fileend, std::string dir, int n)
 
void CopyFile (std::string file1end, std::string file2end)
 
bool CheckEigIsStationary (bool reset=false)
 
bool CheckIfAtNeutralPoint (void)
 
void UpdateAlpha (int n)
 
void UpdateWaveForceMag (int n)
 
void UpdateDAlphaDWaveForceMag (NekDouble alphainit)
 
void AppendEvlToFile (std::string file, int n)
 
void AppendEvlToFile (std::string file, NekDouble WaveForceMag)
 
int GetIterStart ()
 
int GetIterEnd ()
 
VWIIterationType GetVWIIterationType (void)
 
int GetNOuterIterations (void)
 
int GetMaxOuterIterations (void)
 
NekDouble GetAlpha (void)
 
NekDouble GetAlphaStep (void)
 
NekDouble GetWaveForceMag (void)
 
NekDouble GetWaveForceMagStep (void)
 
NekDouble GetDAlphaDWaveForceMag (void)
 
int GetMaxWaveForceMagIter (void)
 
NekDouble GetEigRelTol (void)
 
int GetMinInnerIterations (void)
 
NekDouble GetPrevAlpha (void)
 
void SetAlpha (NekDouble alpha)
 
void SetWaveForceMag (NekDouble mag)
 
void SetEigRelTol (NekDouble tol)
 
void SetAlphaStep (NekDouble step)
 
void SetMinInnerIterations (int niter)
 
void SetPrevAlpha (NekDouble alpha)
 
bool IfIterInterface (void)
 
Array< OneD, int > GetReflectionIndex (void)
 
void FileRelaxation (int reg)
 

Private Attributes

int m_iterStart
 
int m_iterEnd
 
int m_nOuterIterations
 
int m_maxOuterIterations
 
int m_minInnerIterations
 
int m_maxWaveForceMagIter
 
bool m_deltaFcnApprox
 
bool m_useLinfPressureNorm
 
bool m_moveMeshToCriticalLayer
 
NekDouble m_deltaFcnDecay
 
Array< OneD, NekDouble > m_waveForceMag
 
NekDouble m_waveForceMagStep
 
NekDouble m_rollForceScale
 
Array< OneD, NekDouble > m_leading_real_evl
 
Array< OneD, NekDouble > m_leading_imag_evl
 < Leading real eigenvalue More...
 
Array< OneD, NekDouble > m_alpha
 < Leading imaginary eigenvalue More...
 
NekDouble m_alphaStep
 
NekDouble m_neutralPointTol
 
NekDouble m_eigRelTol
 
NekDouble m_vwiRelaxation
 
NekDouble m_dAlphaDWaveForceMag
 
NekDouble m_prevAlpha
 
bool m_iterinterface
 
VWIIterationType m_VWIIterationType
 
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_waveVelocities
 
MultiRegions::ExpListSharedPtr m_wavePressure
 
Array< OneD, Array< OneD,
NekDouble > > 		m_vwiForcing
 
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
 
Array< OneD, Array< OneD,
NekDouble > > 		m_bcsForcing
 
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_streakField
 
Array< OneD,
MultiRegions::ExpListSharedPtr > 		m_rollField
 
std::string m_sessionName
 
LibUtilities::SessionReaderSharedPtr m_sessionVWI
 
LibUtilities::SessionReaderSharedPtr m_sessionRoll
 
EquationSystemSharedPtr m_solverRoll
 
LibUtilities::SessionReaderSharedPtr m_sessionStreak
 
LibUtilities::SessionReaderSharedPtr m_sessionWave
 

Detailed Description

Definition at line 72 of file VortexWaveInteraction.h.

Constructor & Destructor Documentation

Nektar::VortexWaveInteraction::VortexWaveInteraction ( int  argc,
char *  argv[] 
)

Definition at line 47 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::eFixedAlpha, Nektar::eFixedAlphaWaveForcing, Nektar::eFixedWaveForcing, Nektar::eVWIIterationTypeSize, Nektar::SolverUtils::GetEquationSystemFactory(), m_alpha, m_alphaStep, m_dAlphaDWaveForceMag, m_deltaFcnApprox, m_deltaFcnDecay, m_eigRelTol, m_iterEnd, m_iterinterface, m_iterStart, m_leading_imag_evl, m_leading_real_evl, m_maxOuterIterations, m_maxWaveForceMagIter, m_minInnerIterations, m_moveMeshToCriticalLayer, m_neutralPointTol, m_nOuterIterations, m_rollForceScale, m_sessionName, m_sessionRoll, m_sessionStreak, m_sessionVWI, m_sessionWave, m_solverRoll, m_useLinfPressureNorm, m_vwiForcing, m_VWIIterationType, m_vwiRelaxation, m_waveForceMag, m_waveForceMagStep, UpdateAlpha(), and Nektar::VWIIterationTypeMap.

47  :
48  m_nOuterIterations(0)
49  {
50 
51  int storesize;// number of previous iterations to store
52 
53  m_sessionName = argv[argc-1];
54  string meshfile = m_sessionName + ".xml";
55 
56  LibUtilities::SessionReaderSharedPtr session;
57 
58  std::string VWICondFile(argv[argc-1]);
59  VWICondFile += "_VWI.xml";
60  std::vector<std::string> VWIFilenames;
61  VWIFilenames.push_back(meshfile);
62  VWIFilenames.push_back(VWICondFile);
63 
64  // Create Incompressible NavierStokesSolver session reader.
65  m_sessionVWI = LibUtilities::SessionReader::CreateInstance(argc, argv,
66  VWIFilenames);
67 
68  m_sessionVWI->LoadParameter("AlphaStep", m_alphaStep,0.05);
69  m_sessionVWI->LoadParameter("OuterIterationStoreSize",storesize,10);
70 
71  //set a low tol for interfaceVWI
72  m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol,1e-3);
73 
74 
75  m_sessionVWI->LoadParameter("NeutralPointTolerance", m_neutralPointTol,1e-4);
76  m_sessionVWI->LoadParameter("MaxOuterIterations", m_maxOuterIterations,100);
77 
78  m_sessionVWI->LoadParameter("StartIteration",m_iterStart, 0);
79  m_sessionVWI->LoadParameter("EndIteration", m_iterEnd, 0);
80 
81  m_sessionVWI->LoadParameter("WaveForceMagStep", m_waveForceMagStep,0.01);
82  m_sessionVWI->LoadParameter("DAlphaDWaveForceMag", m_dAlphaDWaveForceMag,0.0);
83  m_sessionVWI->LoadParameter("MaxWaveForceMagIter",m_maxWaveForceMagIter,1);
84  m_sessionVWI->LoadParameter("RollForceScale", m_rollForceScale,1.0);
85 
86  if(m_sessionVWI->DefinesSolverInfo("DeltaFcnApprox"))
87  {
88  m_deltaFcnApprox = true;
89  m_sessionVWI->LoadParameter("DeltaFcnDecay", m_deltaFcnDecay,1.0/500);
90  }
91  else
92  {
93  m_deltaFcnApprox = false;
94  m_deltaFcnDecay = 0.0;
95  }
96 
97  if(m_sessionVWI->DefinesSolverInfo("LinfPressureNorm"))
98  {
99  m_useLinfPressureNorm = true;
100  }
101  else
102  {
103  m_useLinfPressureNorm = false;
104  }
105 
106  if( m_sessionVWI->DefinesSolverInfo("INTERFACE") )
107  {
108  m_iterinterface = true;
109  }
110  else
111  {
112  m_iterinterface = false;
113  }
114 
115  if(m_sessionVWI->DefinesSolverInfo("MoveMeshToCriticalLayer"))
116  {
117  m_moveMeshToCriticalLayer = true;
118  }
119  else
120  {
121  m_moveMeshToCriticalLayer = false;
122  }
123 
124  m_alpha = Array<OneD, NekDouble> (storesize);
125  m_alpha[0] = m_sessionVWI->GetParameter("Alpha");
126  m_waveForceMag = Array<OneD, NekDouble> (storesize);
127  m_waveForceMag[0] = m_sessionVWI->GetParameter("WaveForceMag");
128 
129  m_leading_real_evl = Array<OneD, NekDouble> (storesize);
130  m_leading_imag_evl = Array<OneD, NekDouble> (storesize);
131 
132  if(m_sessionVWI->DefinesParameter("Relaxation"))
133  {
134  m_vwiRelaxation = m_sessionVWI->GetParameter("Relaxation");
135  // fix minimum number of iterations to be number of
136  // iterations required to make contribution of innitial
137  // forcing to 0.1
138  m_minInnerIterations = (int) (log(0.1)/log(m_vwiRelaxation));
139  }
140  else
141  {
142  m_vwiRelaxation = 0.0;
143  m_minInnerIterations = 1;
144  }
145 
146  // Initialise NS Roll solver
147  std::string IncCondFile(argv[argc-1]);
148  IncCondFile += "_IncNSCond.xml";
149  std::vector<std::string> IncNSFilenames;
150  IncNSFilenames.push_back(meshfile);
151  IncNSFilenames.push_back(IncCondFile);
152 
153  // Create Incompressible NavierStokesSolver session reader.
154  m_sessionRoll = LibUtilities::SessionReader::CreateInstance(argc, argv, IncNSFilenames,
155  m_sessionVWI->GetComm());
156  std::string vEquation = m_sessionRoll->GetSolverInfo("SolverType");
157  m_solverRoll = GetEquationSystemFactory().CreateInstance(vEquation, m_sessionRoll);
158  m_solverRoll->PrintSummary(cout);
159 
160 
161  if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
162  {
163  m_sessionVWI->LoadParameter("EigenvalueRelativeTol", m_eigRelTol,1e-2);
164  }
165 
166  int ncoeffs = m_solverRoll->UpdateFields()[0]->GetNcoeffs();
167  m_vwiForcing = Array<OneD, Array<OneD, NekDouble> > (4);
168  m_vwiForcing[0] = Array<OneD, NekDouble> (4*ncoeffs);
169  for(int i = 1; i < 4; ++i)
170  {
171  m_vwiForcing[i] = m_vwiForcing[i-1] + ncoeffs;
172  }
173 
174  // Fill forcing into m_vwiForcing
175  // Has forcing even been initialised yet?
176 // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[0]->GetCoeffs(),1,m_vwiForcing[0],1);
177 // Vmath::Vcopy(ncoeffs,m_solverRoll->UpdateForces()[1]->GetCoeffs(),1,m_vwiForcing[1],1);
178 
179 
180  // Create AdvDiff Streak solver
181  std::string AdvDiffCondFile(argv[argc-1]);
182  AdvDiffCondFile += "_AdvDiffCond.xml";
183  std::vector<std::string> AdvDiffFilenames;
184  AdvDiffFilenames.push_back(meshfile);
185  AdvDiffFilenames.push_back(AdvDiffCondFile);
186 
187  // Create AdvDiffusion session reader.
188  m_sessionStreak = LibUtilities::SessionReader::CreateInstance(argc, argv, AdvDiffFilenames, m_sessionVWI->GetComm());
189 
190  // Initialise LinNS solver
191  std::string LinNSCondFile(argv[argc-1]);
192  LinNSCondFile += "_LinNSCond.xml";
193  std::vector<std::string> LinNSFilenames;
194  LinNSFilenames.push_back(meshfile);
195  LinNSFilenames.push_back(LinNSCondFile);
196 
197  // Create Linearised NS stability session reader.
198  m_sessionWave = LibUtilities::SessionReader::CreateInstance(argc, argv, LinNSFilenames, m_sessionVWI->GetComm());
199 
200  // Set the initial beta value in stability to be equal to VWI file
201  std::string LZstr("LZ");
202  NekDouble LZ = 2*M_PI/m_alpha[0];
203  cout << "Setting LZ in Linearised solver to " << LZ << endl;
204  m_sessionWave->SetParameter(LZstr,LZ);
205 
206  // Check for iteration type
207  if(m_sessionVWI->DefinesSolverInfo("VWIIterationType"))
208  {
209  std::string IterationTypeStr = m_sessionVWI->GetSolverInfo("VWIIterationType");
210  for(int i = 0; i < (int) eVWIIterationTypeSize; ++i)
211  {
212  if(m_solverRoll->NoCaseStringCompare(VWIIterationTypeMap[i],IterationTypeStr) == 0 )
213  {
214  m_VWIIterationType = (VWIIterationType)i;
215  break;
216  }
217  }
218  }
219  else
220  {
221  m_VWIIterationType = eFixedAlphaWaveForcing;
222  }
223 
224 
225 
226  // Check for restart
227  bool restart;
228  m_sessionVWI->MatchSolverInfo("RestartIteration","True",restart,false);
229  if(restart)
230  {
231  switch(m_VWIIterationType)
232  {
233  case eFixedAlpha:
234  break;
235  case eFixedWaveForcing:
236  {
237  FILE *fp;
238  // Check for OuterIter.his file to read
239  if((fp = fopen("OuterIter.his","r")))
240  {
241  char buf[BUFSIZ];
242  std::vector<NekDouble> Alpha, Growth, Phase;
243  NekDouble alpha,growth,phase;
244  while(fgets(buf,BUFSIZ,fp))
245  {
246  sscanf(buf,"%*d:%lf%lf%lf",&alpha,&growth,&phase);
247  Alpha.push_back(alpha);
248  Growth.push_back(growth);
249  Phase.push_back(phase);
250  }
251 
252  m_nOuterIterations = Alpha.size();
253 
254  int nvals = std::min(m_nOuterIterations,(int)m_alpha.num_elements());
255 
256  for(int i = 0; i < nvals; ++i)
257  {
258  m_alpha[nvals-1-i] = Alpha[m_nOuterIterations-nvals+i];
259  m_leading_real_evl[nvals-1-i] = Growth[m_nOuterIterations-nvals+i];
260  m_leading_imag_evl[nvals-1-i] = Phase [m_nOuterIterations-nvals+i];
261  }
262 
263  UpdateAlpha(m_nOuterIterations++);
264  }
265  else
266  {
267  cout << " No File OuterIter.his to restart from" << endl;
268  }
269  }
270  break;
271  case eFixedAlphaWaveForcing:
272  {
273  string nstr = boost::lexical_cast<std::string>(m_iterStart);
274  cout << "Restarting from iteration " << m_iterStart << endl;
275  std::string rstfile = "cp -f Save/" + m_sessionName + ".rst." + nstr + " " + m_sessionName + ".rst";
276  cout << " " << rstfile << endl;
277  if(system(rstfile.c_str()))
278  {
279  ASSERTL0(false,rstfile.c_str());
280  }
281  std::string vwifile = "cp -f Save/" + m_sessionName + ".vwi." + nstr + " " + m_sessionName + ".vwi";
282  cout << " " << vwifile << endl;
283  if(system(vwifile.c_str()))
284  {
285  ASSERTL0(false,vwifile.c_str());
286  }
287  }
288  break;
289  default:
290  ASSERTL0(false,"Unknown VWIITerationType in restart");
291  }
292  }
293 
294  // Check for ConveredSoln to update DAlphaDWaveForce
295  {
296  FILE *fp;
297  if((fp = fopen("ConvergedSolns","r")))
298  {
299  char buf[BUFSIZ];
300  std::vector<NekDouble> WaveForce, Alpha;
301  NekDouble waveforce,alpha;
302  while(fgets(buf,BUFSIZ,fp))
303  {
304  sscanf(buf,"%*d:%lf%lf",&waveforce,&alpha);
305  WaveForce.push_back(waveforce);
306  Alpha.push_back(alpha);
307  }
308 
309  if(Alpha.size() > 1)
310  {
311  int min_i = 0;
312  NekDouble min_alph = fabs(m_alpha[0]-Alpha[min_i]);
313  // find nearest point
314  for(int i = 1; i < Alpha.size(); ++i)
315  {
316  if(fabs(m_alpha[0]-Alpha[min_i]) < min_alph)
317  {
318  min_i = i;
319  min_alph = fabs(m_alpha[0]-Alpha[min_i]);
320  }
321  }
322 
323  // find next nearest point
324  int min_j = (min_i == 0)? 1:0;
325  min_alph = fabs(m_alpha[0]-Alpha[min_j]);
326  for(int i = 0; i < Alpha.size(); ++i)
327  {
328  if(i != min_i)
329  {
330  if(fabs(m_alpha[0]-Alpha[min_j]) < min_alph)
331  {
332  min_j = i;
333  min_alph = fabs(m_alpha[0]-Alpha[min_j]);
334  }
335  }
336  }
337 
338  if(fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
339  {
340  m_dAlphaDWaveForceMag = (Alpha[min_i]-Alpha[min_j])/(WaveForce[min_i]-WaveForce[min_j]);
341  }
342  }
343  }
344  }
345  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::VortexWaveInteraction::m_sessionRoll
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Definition: VortexWaveInteraction.h:272
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::VortexWaveInteraction::m_solverRoll
EquationSystemSharedPtr m_solverRoll
Definition: VortexWaveInteraction.h:273
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:241
Nektar::LibUtilities::SessionReader::CreateInstance
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
Definition: SessionReader.h:140
Nektar::VortexWaveInteraction::m_sessionWave
LibUtilities::SessionReaderSharedPtr m_sessionWave
Definition: VortexWaveInteraction.h:276
Nektar::VWIIterationTypeMap
const std::string VWIIterationTypeMap[]
Definition: VortexWaveInteraction.h:64
Nektar::VortexWaveInteraction::m_waveForceMag
Array< OneD, NekDouble > m_waveForceMag
Definition: VortexWaveInteraction.h:236
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:83
Nektar::VortexWaveInteraction::m_vwiForcing
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
Definition: VortexWaveInteraction.h:260
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
Nektar::VortexWaveInteraction::m_sessionStreak
LibUtilities::SessionReaderSharedPtr m_sessionStreak
Definition: VortexWaveInteraction.h:275
Nektar::VortexWaveInteraction::~VortexWaveInteraction ( void  )

Definition at line 348 of file VortexWaveInteraction.cpp.

References m_sessionVWI.

349  {
350  m_sessionVWI->Finalise();
351  }
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270

Member Function Documentation

void Nektar::VortexWaveInteraction::AppendEvlToFile ( std::string  file,
int  n 
)

Definition at line 1062 of file VortexWaveInteraction.cpp.

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

Referenced by DoFixedForcingIteration(), and main().

1063  {
1064  FILE *fp;
1065  fp = fopen(file.c_str(),"a");
1066  fprintf(fp, "%d: %lf %16.12le %16.12le\n",n, m_alpha[0], m_leading_real_evl[0],m_leading_imag_evl[0]);
1067  fclose(fp);
1068  }
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:241
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
void Nektar::VortexWaveInteraction::AppendEvlToFile ( std::string  file,
NekDouble  WaveForceMag 
)

Definition at line 1070 of file VortexWaveInteraction.cpp.

References m_alpha, m_leading_imag_evl, and m_leading_real_evl.

1071  {
1072  FILE *fp;
1073  fp = fopen(file.c_str(),"a");
1074  fprintf(fp, "%lf %lf %16.12le %16.12le\n",WaveForceMag, m_alpha[0], m_leading_real_evl[0],m_leading_imag_evl[0]);
1075  fclose(fp);
1076  }
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:241
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
void Nektar::VortexWaveInteraction::CalcL2ToLinfPressure ( void  )

Definition at line 965 of file VortexWaveInteraction.cpp.

References ExecuteWave(), Vmath::Fill(), m_wavePressure, m_waveVelocities, npts, Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), and Vmath::Vvtvp().

Referenced by main().

966  {
967 
968  ExecuteWave();
969 
970  m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),
971  m_wavePressure->GetPlane(0)->UpdatePhys());
972  m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),
973  m_wavePressure->GetPlane(1)->UpdatePhys());
974 
975  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
976  NekDouble Linf;
977  Array<OneD, NekDouble> val(2*npts,0.0);
978 
979  Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1,val,1);
980  Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1,val,1,val,1);
981  cout << "int P^2 " << m_wavePressure->GetPlane(0)->PhysIntegral(val) << endl;
982  Vmath::Vsqrt(npts,val,1,val,1);
983 
984 
985  Linf = Vmath::Vmax(npts,val,1);
986  cout << "Linf: " << Linf << endl;
987 
988  NekDouble l2,norm;
989  l2 = m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
990  norm = l2*l2;
991  l2 = m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
992  norm += l2*l2;
993 
994 
995  Vmath::Fill(npts,1.0,val,1);
996  NekDouble area = m_waveVelocities[0]->GetPlane(0)->PhysIntegral(val);
997 
998  l2 = sqrt(norm/area);
999 
1000  cout << "L2: " << l2 << endl;
1001 
1002  cout << "Ratio Linf/L2: "<< Linf/l2 << endl;
1003  }
Nektar::VortexWaveInteraction::m_wavePressure
MultiRegions::ExpListSharedPtr m_wavePressure
Definition: VortexWaveInteraction.h:258
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Nektar::VortexWaveInteraction::m_waveVelocities
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
Definition: VortexWaveInteraction.h:257
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
void Nektar::VortexWaveInteraction::CalcNonLinearWaveForce ( void  )

Definition at line 560 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::SpatialDomains::eNeumann, Vmath::Fill(), GetReflectionIndex(), m_alpha, m_deltaFcnApprox, m_deltaFcnDecay, m_sessionName, m_sessionRoll, m_sessionVWI, m_solverRoll, m_streakField, m_useLinfPressureNorm, m_vwiForcing, m_vwiRelaxation, m_waveForceMag, m_wavePressure, m_waveVelocities, Vmath::Neg(), npts, Vmath::Sadd(), Vmath::Smul(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmax(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

Referenced by DoFixedForcingIteration(), ExecuteLoop(), and main().

561  {
562  if(m_sessionRoll->DefinesSolverInfo("INTERFACE") )
563  {
564  static int cnt=0;
565  string wavefile = m_sessionName +".fld";
566  string movedmesh = m_sessionName + "_advPost_moved.xml";
567  string filestreak = m_sessionName + "_streak.fld";
568  char c[16]="";
569  sprintf(c,"%d",cnt);
570  char c_alpha[16]="";
571  sprintf(c_alpha,"%f",m_alpha[0]);
572  string syscall;
573  if( m_sessionVWI->GetSolverInfo("INTERFACE")=="phase" )
574  {
575  string filePost = m_sessionName + "_advPost.xml";
576  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs "
577  + filePost +" "+
578  "meshhalf_pos_Spen_stability_moved.fld meshhalf_pos_Spen_advPost_moved.fld "
579  +c_alpha +" > data_alpha0";
580  cout<<syscall.c_str()<<endl;
581  if(system(syscall.c_str()))
582  {
583  ASSERTL0(false,syscall.c_str());
584  }
585 
586  syscall = "cp -f meshhalf_pos_Spen_stability_moved_u_5.bc "+m_sessionName+"_u_5.bc";
587  cout<<syscall.c_str()<<endl;
588  if(system(syscall.c_str()))
589  {
590  ASSERTL0(false,syscall.c_str());
591  }
592  syscall = "cp -f meshhalf_pos_Spen_stability_moved_v_5.bc "+m_sessionName+"_v_5.bc";
593  cout<<syscall.c_str()<<endl;
594  if(system(syscall.c_str()))
595  {
596  ASSERTL0(false,syscall.c_str());
597  }
598  }
599  else
600  {
601  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs "
602  + movedmesh + " " + wavefile + " " + filestreak + " "+c_alpha +" > datasub_"+c;
603  cout<<syscall.c_str()<<endl;
604  if(system(syscall.c_str()))
605  {
606  ASSERTL0(false,syscall.c_str());
607  }
608  }
609 
610 
611 
612  //relaxation for different alpha values? does it make sense?
613 
614  //save the wave
615  string wave_subalp = m_sessionName + "_wave_subalp_"+c+".fld";
616  syscall = "cp -f " + wavefile + " " + wave_subalp;
617  cout<<syscall.c_str()<<endl;
618  if(system(syscall.c_str()))
619  {
620  ASSERTL0(false,syscall.c_str());
621  }
622  //FileRelaxation(3);
623  cnt++;
624 
625  }
626  else
627  {
628  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
629  int ncoeffs = m_waveVelocities[0]->GetPlane(0)->GetNcoeffs();
630  Array<OneD, NekDouble> val(npts), der1(2*npts);
631  Array<OneD, NekDouble> der2 = der1 + npts;
632  Array<OneD, NekDouble> streak;
633  static int projectfield = -1;
634 
635  if(m_deltaFcnApprox)
636  {
637  streak = Array<OneD, NekDouble> (npts);
638  m_streakField[0]->BwdTrans(m_streakField[0]->GetCoeffs(), streak);
639  }
640 
641  // Set project field to be first field that has a Neumann
642  // boundary since this not impose any condition on the vertical boundaries
643  // Othersise set to zero.
644  if(projectfield == -1)
645  {
646  Array<OneD, const SpatialDomains::BoundaryConditionShPtr > BndConds;
647 
648  for(int i = 0; i < m_waveVelocities.num_elements(); ++i)
649  {
650  BndConds = m_waveVelocities[i]->GetBndConditions();
651  for(int j = 0; j < BndConds.num_elements(); ++j)
652  {
653  if(BndConds[j]->GetBoundaryConditionType() == SpatialDomains::eNeumann)
654  {
655  projectfield = i;
656  break;
657  }
658  }
659  if(projectfield != -1)
660  {
661  break;
662  }
663  }
664  if(projectfield == -1)
665  {
666  cout << "using first field to project non-linear forcing which imposes a Dirichlet condition" << endl;
667  projectfield = 0;
668  }
669  }
670 
671  // Shift m_vwiForcing in case of relaxation
672  Vmath::Vcopy(ncoeffs,m_vwiForcing[0],1,m_vwiForcing[2],1);
673  Vmath::Vcopy(ncoeffs,m_vwiForcing[1],1,m_vwiForcing[3],1);
674 
675  // determine inverse of area normalised field.
676  m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(), m_wavePressure->GetPlane(0)->UpdatePhys());
677  m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(), m_wavePressure->GetPlane(1)->UpdatePhys());
678  NekDouble invnorm;
679 
680  if(m_useLinfPressureNorm)
681  {
682  Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1,der1,1);
683  Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1,der1,1,der1,1);
684  Vmath::Vsqrt(npts,der1,1,der1,1);
685 
686  NekDouble Linf = Vmath::Vmax(npts,der1,1);
687 
688  invnorm = 1.0/Linf;
689  }
690  else
691  {
692  // Determine normalisation of pressure so that |P|/A = 1
693  NekDouble l2;
694  l2 = m_wavePressure->GetPlane(0)->L2(m_wavePressure->GetPlane(0)->GetPhys());
695  invnorm = l2*l2;
696  l2 = m_wavePressure->GetPlane(1)->L2(m_wavePressure->GetPlane(1)->GetPhys());
697  invnorm += l2*l2;
698  Vmath::Fill(2*npts,1.0,der1,1);
699  NekDouble area = m_waveVelocities[0]->GetPlane(0)->PhysIntegral(der1);
700  cout << "Area: " << area << endl;
701  invnorm = sqrt(area/invnorm);
702  }
703 
704  // Get hold of arrays.
705  m_waveVelocities[0]->GetPlane(0)->BwdTrans(m_waveVelocities[0]->GetPlane(0)->GetCoeffs(),m_waveVelocities[0]->GetPlane(0)->UpdatePhys());
706  Array<OneD, NekDouble> u_real = m_waveVelocities[0]->GetPlane(0)->UpdatePhys();
707  Vmath::Smul(npts,invnorm,u_real,1,u_real,1);
708  m_waveVelocities[0]->GetPlane(1)->BwdTrans(m_waveVelocities[0]->GetPlane(1)->GetCoeffs(),m_waveVelocities[0]->GetPlane(1)->UpdatePhys());
709  Array<OneD, NekDouble> u_imag = m_waveVelocities[0]->GetPlane(1)->UpdatePhys();
710  Vmath::Smul(npts,invnorm,u_imag,1,u_imag,1);
711  m_waveVelocities[1]->GetPlane(0)->BwdTrans(m_waveVelocities[1]->GetPlane(0)->GetCoeffs(),m_waveVelocities[1]->GetPlane(0)->UpdatePhys());
712  Array<OneD, NekDouble> v_real = m_waveVelocities[1]->GetPlane(0)->UpdatePhys();
713  Vmath::Smul(npts,invnorm,v_real,1,v_real,1);
714  m_waveVelocities[1]->GetPlane(1)->BwdTrans(m_waveVelocities[1]->GetPlane(1)->GetCoeffs(),m_waveVelocities[1]->GetPlane(1)->UpdatePhys());
715  Array<OneD, NekDouble> v_imag = m_waveVelocities[1]->GetPlane(1)->UpdatePhys();
716  Vmath::Smul(npts,invnorm,v_imag,1,v_imag,1);
717 
718  // normalise wave for output
719  Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[0]->UpdateCoeffs(),1,m_waveVelocities[0]->UpdateCoeffs(),1);
720  Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[1]->UpdateCoeffs(),1,m_waveVelocities[1]->UpdateCoeffs(),1);
721  Vmath::Smul(2*ncoeffs,invnorm,m_waveVelocities[2]->UpdateCoeffs(),1,m_waveVelocities[2]->UpdateCoeffs(),1);
722 
723  // dump field
724  {
725  std::vector<std::string> variables(3);
726  std::vector<Array<OneD, NekDouble> > outfield(3);
727  variables[0] = "u_w";
728  variables[1] = "v_w";
729  variables[2] = "w_w";
730  outfield[0] = m_waveVelocities[0]->UpdateCoeffs();
731  outfield[1] = m_waveVelocities[1]->UpdateCoeffs();
732  outfield[2] = m_waveVelocities[2]->UpdateCoeffs();
733  std::string outname = m_sessionName + "_wave.fld";
734  m_solverRoll->WriteFld(outname, m_waveVelocities[0], outfield, variables);
735  }
736 
737 #if 1
738  int ncoeffs_p = m_wavePressure->GetPlane(0)->GetNcoeffs();
739  Vmath::Smul(ncoeffs_p,invnorm,m_wavePressure->GetPlane(0)->UpdateCoeffs(),1,m_wavePressure->GetPlane(0)->UpdateCoeffs(),1);
740  Vmath::Smul(ncoeffs_p,invnorm,m_wavePressure->GetPlane(1)->UpdateCoeffs(),1,m_wavePressure->GetPlane(1)->UpdateCoeffs(),1);
741 #else
742  m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),m_wavePressure->GetPlane(0)->UpdatePhys());
743  Vmath::Smul(npts,invnorm,m_wavePressure->GetPlane(0)->UpdatePhys(),1,m_wavePressure->GetPlane(0)->UpdatePhys(),1);
744  m_wavePressure->GetPlane(0)->FwdTrans(m_wavePressure->GetPlane(0)->UpdatePhys(),m_wavePressure->GetPlane(0)->UpdateCoeffs());
745 
746  m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),m_wavePressure->GetPlane(1)->UpdatePhys());
747  Vmath::Smul(npts,invnorm,m_wavePressure->GetPlane(1)->UpdatePhys(),1,m_wavePressure->GetPlane(1)->UpdatePhys(),1);
748  m_wavePressure->GetPlane(1)->FwdTrans(m_wavePressure->GetPlane(1)->UpdatePhys(),m_wavePressure->GetPlane(1)->UpdateCoeffs());
749 #endif
750 
751  // Calculate non-linear terms for x and y directions
752  // d/dx(u u* + u* u)
753  Vmath::Vmul (npts,u_real,1,u_real,1,val,1);
754  Vmath::Vvtvp(npts,u_imag,1,u_imag,1,val,1,val,1);
755  Vmath::Smul (npts,2.0,val,1,val,1);
756  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
757 
758  // d/dy(v u* + v* u)
759  Vmath::Vmul (npts,u_real,1,v_real,1,val,1);
760  Vmath::Vvtvp(npts,u_imag,1,v_imag,1,val,1,val,1);
761  Vmath::Smul (npts,2.0,val,1,val,1);
762  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
763 
764  Vmath::Vadd(npts,der1,1,der2,1,der1,1);
765 #if 1
766  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[0]);
767 #else
768  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,m_vwiForcing[0]);
769 #endif
770  Vmath::Smul(ncoeffs,-m_waveForceMag[0],m_vwiForcing[0],1,m_vwiForcing[0],1);
771 
772  // d/dx(u v* + u* v)
773  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(0,val,der1);
774 
775  // d/dy(v v* + v* v)
776  Vmath::Vmul(npts,v_real,1,v_real,1,val,1);
777  Vmath::Vvtvp(npts,v_imag,1,v_imag,1,val,1,val,1);
778  Vmath::Smul (npts,2.0,val,1,val,1);
779  m_waveVelocities[0]->GetPlane(0)->PhysDeriv(1,val,der2);
780 
781  Vmath::Vadd(npts,der1,1,der2,1,der1,1);
782 
783 #if 1
784  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[1]);
785 #else
786  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der1,m_vwiForcing[1]);
787 #endif
788 
789  Vmath::Smul(ncoeffs,-m_waveForceMag[0],m_vwiForcing[1],1,m_vwiForcing[1],1);
790 
791  //by default the symmetrization is on
792  bool symm=true;
793  m_sessionVWI->MatchSolverInfo("Symmetrization","True",symm,true);
794 #if 0
795  if(symm== true )
796  {
797 
798  // Symmetrise forcing
799  //-> Get coordinates
800  Array<OneD, NekDouble> coord(2);
801  Array<OneD, NekDouble> coord_x(npts);
802  Array<OneD, NekDouble> coord_y(npts);
803 
804  //-> Impose symmetry (x -> -x + Lx/2, y-> -y) on coordinates
805  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
806  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
807  Vmath::Neg(npts,coord_x,1);
808  Vmath::Sadd(npts,xmax,coord_x,1,coord_x,1);
809  Vmath::Neg(npts,coord_y,1);
810 
811  int i, physoffset;
812 
813  //-> Obtain list of expansion element ids for each point.
814  Array<OneD, int> Eid(npts);
815  // This search may not be necessary every iteration
816  for(i = 0; i < npts; ++i)
817  {
818  coord[0] = coord_x[i];
819  coord[1] = coord_y[i];
820 
821  // Note this will not quite be symmetric.
822  Eid[i] = m_waveVelocities[0]->GetPlane(0)->GetExpIndex(coord,1e-6);
823  }
824 
825  // Interpolate field 0
826  m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[0],der1);
827  for(i = 0; i < npts; ++i)
828  {
829  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
830  coord[0] = coord_x[i];
831  coord[1] = coord_y[i];
832  der2 [i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
833  }
834  //-> Average field 0
835  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
836  Vmath::Smul(npts,0.5,der2,1,der2,1);
837 #if 1
838  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[0]);
839 #else
840  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForcing[0]);
841 #endif
842 
843  //-> Interpoloate field 1
844  m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[1],der1);
845  for(i = 0; i < npts; ++i)
846  {
847  physoffset = m_waveVelocities[0]->GetPlane(0)->GetPhys_Offset(Eid[i]);
848  coord[0] = coord_x[i];
849  coord[1] = coord_y[i];
850  der2[i] = m_waveVelocities[0]->GetPlane(0)->GetExp(Eid[i])->PhysEvaluate(coord, der1 + physoffset);
851  }
852 
853  //-> Average field 1
854  Vmath::Vsub(npts,der1,1,der2,1,der2,1);
855  Vmath::Smul(npts,0.5,der2,1,der2,1);
856 #if 1
857  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der2,m_vwiForcing[1]);
858 #else
859  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
860 #endif
861  }
862 #else
863  int i;
864  if(symm== true )
865  {
866  cout<<"symmetrization is active"<<endl;
867  static Array<OneD, int> index = GetReflectionIndex();
868 
869  m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[0],der1);
870  for(i = 0; i < npts; ++i)
871  {
872  if(index[i] != -1)
873  {
874  val[i] = 0.5*(der1[i] - der1[index[i]]);
875  }
876  }
877 #if 1
878  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(val,m_vwiForcing[0]);
879 #else
880  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, m_vwiForcing[0]);
881 #endif
882 
883  m_waveVelocities[0]->GetPlane(0)->BwdTrans_IterPerExp(m_vwiForcing[1],der2);
884  for(i = 0; i < npts; ++i)
885  {
886  if(index[i] != -1)
887  {
888  val[i] = 0.5*(der2[i] - der2[index[i]]);
889  }
890  }
891 #if 1
892  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(val,m_vwiForcing[1]);
893 #else
894  m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(val, m_vwiForcing[1]);
895 #endif
896  }
897 
898 
899  Vmath::Vmul(npts,der1,1,der1,1,val,1);
900  Vmath::Vvtvp(npts,der2,1,der2,1,val,1,val,1);
901  Vmath::Vsqrt(npts,val,1,val,1);
902  cout << "F_Linf: " << Vmath::Vmax(npts,val,1) << endl;
903 
904 #endif
905 
906  if(m_vwiRelaxation)
907  {
908  Vmath::Smul(ncoeffs,1.0-m_vwiRelaxation,
909  m_vwiForcing[0],1,m_vwiForcing[0],1);
910  Vmath::Svtvp(ncoeffs,m_vwiRelaxation,m_vwiForcing[2],1,
911  m_vwiForcing[0],1,m_vwiForcing[0],1);
912 
913  Vmath::Smul(ncoeffs,1.0-m_vwiRelaxation,
914  m_vwiForcing[1],1,m_vwiForcing[1],1);
915  Vmath::Svtvp(ncoeffs,m_vwiRelaxation,m_vwiForcing[3],1,
916  m_vwiForcing[1],1,m_vwiForcing[1],1);
917  }
918 
919  if(m_deltaFcnApprox)
920  {
921  for(int j = 0; j < 2; ++j)
922  {
923 
924  m_waveVelocities[projectfield]->GetPlane(0)->BwdTrans(m_vwiForcing[j],der1);
925  for(int i = 0; i < npts; ++i)
926  {
927  der1[i] *= exp(-streak[i]*streak[i]/m_deltaFcnDecay);
928  }
929  m_waveVelocities[projectfield]->GetPlane(0)->FwdTrans(der1,m_vwiForcing[j]);
930  }
931  }
932 
933 
934  // dump output
935  std::vector<std::string> variables(4);
936  std::vector<Array<OneD, NekDouble> > outfield(4);
937  variables[0] = "u"; variables[1] = "v";
938  variables[2] = "pr"; variables[3] = "pi";
939  outfield[0] = m_vwiForcing[0];
940  outfield[1] = m_vwiForcing[1];
941  Array<OneD,NekDouble> soln(npts,0.0);
942  m_wavePressure->GetPlane(0)->BwdTrans(m_wavePressure->GetPlane(0)->GetCoeffs(),m_wavePressure->GetPlane(0)->UpdatePhys());
943  outfield[2] = Array<OneD,NekDouble>(ncoeffs);
944  m_waveVelocities[0]->GetPlane(0)->FwdTrans_IterPerExp(m_wavePressure->GetPlane(0)->GetPhys(),outfield[2]);
945  m_wavePressure->GetPlane(1)->BwdTrans(m_wavePressure->GetPlane(1)->GetCoeffs(),m_wavePressure->GetPlane(1)->UpdatePhys());
946 
947  Vmath::Vmul(npts,m_wavePressure->GetPlane(0)->UpdatePhys(),1,
948  m_wavePressure->GetPlane(0)->UpdatePhys(),1,val,1);
949  Vmath::Vvtvp(npts,m_wavePressure->GetPlane(1)->UpdatePhys(),1,
950  m_wavePressure->GetPlane(1)->UpdatePhys(),1,val,1,val,1);
951  cout << "int P^2: " << m_wavePressure->GetPlane(0)->PhysIntegral(val) << endl;
952  Vmath::Vsqrt(npts,val,1,val,1);
953  cout << "PLinf: " << Vmath::Vmax(npts,val,1) << endl;
954 
955  outfield[3] = Array<OneD,NekDouble>(ncoeffs);
956  m_waveVelocities[1]->GetPlane(0)->FwdTrans_IterPerExp(m_wavePressure->GetPlane(1)->GetPhys(),outfield[3]);
957 
958  std::string outname = m_sessionName + ".vwi";
959 
960  m_solverRoll->WriteFld(outname, m_waveVelocities[0]->GetPlane(0), outfield, variables);
961 
962  }
963  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::VortexWaveInteraction::m_sessionRoll
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Definition: VortexWaveInteraction.h:272
Nektar::VortexWaveInteraction::m_wavePressure
MultiRegions::ExpListSharedPtr m_wavePressure
Definition: VortexWaveInteraction.h:258
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:394
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Nektar::VortexWaveInteraction::m_waveVelocities
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
Definition: VortexWaveInteraction.h:257
Nektar::VortexWaveInteraction::m_solverRoll
EquationSystemSharedPtr m_solverRoll
Definition: VortexWaveInteraction.h:273
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::VortexWaveInteraction::m_waveForceMag
Array< OneD, NekDouble > m_waveForceMag
Definition: VortexWaveInteraction.h:236
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
Nektar::VortexWaveInteraction::m_vwiForcing
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
Definition: VortexWaveInteraction.h:260
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
Nektar::VortexWaveInteraction::m_streakField
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Definition: VortexWaveInteraction.h:265
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
Nektar::VortexWaveInteraction::GetReflectionIndex
Array< OneD, int > GetReflectionIndex(void)
Definition: VortexWaveInteraction.cpp:1925
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
bool Nektar::VortexWaveInteraction::CheckEigIsStationary ( bool  reset = false)

Definition at line 1605 of file VortexWaveInteraction.cpp.

References m_eigRelTol, m_leading_imag_evl, m_leading_real_evl, and m_minInnerIterations.

Referenced by DoFixedForcingIteration().

1606  {
1607  static NekDouble previous_real_evl = -1.0;
1608  static NekDouble previous_imag_evl = -1.0;
1609  static int min_iter = 0;
1610 
1611  if(reset)
1612  {
1613  previous_real_evl = -1.0;
1614  min_iter = 0;
1615  }
1616 
1617  if(previous_real_evl == -1.0 || min_iter < m_minInnerIterations)
1618  {
1619  previous_real_evl = m_leading_real_evl[0];
1620  previous_imag_evl = m_leading_imag_evl[0];
1621  min_iter++;
1622  return false;
1623  }
1624 
1625  cout << "Growth tolerance: " << fabs((m_leading_real_evl[0] - previous_real_evl)/m_leading_real_evl[0]) << endl;
1626  cout << "Phase tolerance: " << fabs((m_leading_imag_evl[0] - previous_imag_evl)/m_leading_imag_evl[0]) << endl;
1627 
1628  // See if real and imaginary growth have converged to with m_eigRelTol
1629  if((fabs((m_leading_real_evl[0] - previous_real_evl)/m_leading_real_evl[0]) < m_eigRelTol)||(fabs(m_leading_real_evl[0]) < 1e-6))
1630  {
1631  previous_real_evl = m_leading_real_evl[0];
1632  previous_imag_evl = m_leading_imag_evl[0];
1633  if((fabs((m_leading_imag_evl[0] - previous_imag_evl)/m_leading_imag_evl[0]) < m_eigRelTol)||(fabs(m_leading_imag_evl[0]) < 1e-6))
1634  {
1635  return true;
1636  }
1637  }
1638  else
1639  {
1640  if(fabs(m_leading_imag_evl[0]) > 1e-2)
1641  {
1642  cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1643  }
1644  previous_real_evl = m_leading_real_evl[0];
1645  previous_imag_evl = m_leading_imag_evl[0];
1646  return false;
1647  }
1648  return false;
1649  }
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:241
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
bool Nektar::VortexWaveInteraction::CheckIfAtNeutralPoint ( void  )

Definition at line 1653 of file VortexWaveInteraction.cpp.

References m_leading_imag_evl, m_leading_real_evl, m_neutralPointTol, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration().

1654  {
1655  bool returnval = false;
1656  if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1657  {
1658  if(m_sessionVWI->GetSolverInfo("INTERFACE")=="phase")
1659  {
1660  if( abs(m_leading_real_evl[0]) < 1e-4 )
1661  {
1662  returnval = true;
1663  }
1664  }
1665  else
1666  {
1667  if( abs(m_leading_real_evl[0]) < 1e-4 && abs(m_leading_imag_evl[0]) <2e-6 )
1668  {
1669  returnval = true;
1670  }
1671  }
1672 
1673  }
1674  else
1675  {
1676  if((m_leading_real_evl[0]*m_leading_real_evl[0] +
1677  m_leading_imag_evl[0]*m_leading_imag_evl[0]) <
1678  m_neutralPointTol*m_neutralPointTol)
1679  {
1680  returnval = true;
1681  }
1682  }
1683  if(fabs(m_leading_imag_evl[0]) > 1e-2)
1684  {
1685  cout << "Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1686  }
1687 
1688  return returnval;
1689  }
Nektar::VortexWaveInteraction::m_sessionRoll
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Definition: VortexWaveInteraction.h:272
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:241
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
void Nektar::VortexWaveInteraction::CopyFile ( std::string  file1end,
std::string  file2end 
)

Definition at line 1050 of file VortexWaveInteraction.cpp.

References ASSERTL0, and m_sessionName.

Referenced by ExecuteRoll(), ExecuteStreak(), and ExecuteWave().

1051  {
1052  string cpfile1 = m_sessionName + file1end;
1053  string cpfile2 = m_sessionName + file2end;
1054  string syscall = "cp -f " + cpfile1 + " " + cpfile2;
1055 
1056  if(system(syscall.c_str()))
1057  {
1058  ASSERTL0(false,syscall.c_str());
1059  }
1060  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
void Nektar::VortexWaveInteraction::ExecuteLoop ( bool  CalcWaveForce = true)

Definition at line 1098 of file VortexWaveInteraction.cpp.

References ASSERTL0, CalcNonLinearWaveForce(), Nektar::LibUtilities::NekManager< KeyType, ValueT, opLessCreator >::DisableManagement(), Nektar::eFixedWaveForcingWithSubIterationOnAlpha, Nektar::LibUtilities::NekManager< KeyType, ValueT, opLessCreator >::EnableManagement(), ExecuteRoll(), ExecuteStreak(), ExecuteWave(), GetVWIIterationType(), m_alpha, m_leading_imag_evl, m_moveMeshToCriticalLayer, m_sessionName, m_sessionRoll, and m_sessionVWI.

Referenced by DoFixedForcingIteration(), and main().

1099  {
1100 
1101 
1102  //the global info has to be written in the
1103  //roll session file to use the interface loop
1104  if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
1105  {
1106  static int cnt=0;
1107  bool skiprollstreak=false;
1108  if(cnt==0 && m_sessionVWI->GetParameter("rollstreakfromit")==1)
1109  {
1110  skiprollstreak =true;
1111  cout<<"skip roll-streak at the first iteration"<<endl;
1112  }
1113 
1114 
1115  if(skiprollstreak != true)
1116  {
1117 
1118  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1119  ExecuteRoll();
1120  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1121 #ifndef _WIN32
1122  sleep(3);
1123 #endif
1124  ExecuteStreak();
1125 #ifndef _WIN32
1126  sleep(3);
1127 #endif
1128  }
1129 
1130  string syscall;
1131  char c[16]="";
1132  string movedmesh = m_sessionName + "_advPost_moved.xml";
1133  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1134  //rewrite the Rollsessionfile (we start from the waleffe forcing)
1135  //string meshbndjumps = m_sessionName +"_bndjumps.xml";
1136  //if(cnt==0)
1137  //{
1138  //take the conditions tag from meshbndjumps and copy into
1139  // the rolls session file
1140  //}
1141 
1142 
1143  sprintf(c,"%d",cnt);
1144  //save old roll solution
1145  string oldroll = m_sessionName +"_roll_"+c +".fld";
1146  syscall = "cp -f " + m_sessionName+"-Base.fld" + " " + oldroll;
1147  cout<<syscall.c_str()<<endl;
1148  if(system(syscall.c_str()))
1149  {
1150  ASSERTL0(false,syscall.c_str());
1151  }
1152  //define file names
1153  string filePost = m_sessionName + "_advPost.xml";
1154  string filestreak = m_sessionName + "_streak.fld";
1155  string filewave = m_sessionName + "_wave.fld";
1156  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1157  string fileinterp = m_sessionName + "_interp.xml";
1158  string interpstreak = m_sessionName +"_interpstreak_"+ c +".fld";
1159  string interwavepressure = m_sessionName +"_wave_p_split_interp_"+ c +".fld";
1160  char alpchar[16]="";
1161 cout<<"alpha = "<<m_alpha[0]<<endl;
1162  sprintf(alpchar, "%f", m_alpha[0]);
1163 
1164 
1165  if( m_sessionVWI->GetSolverInfo("INTERFACE")!="phase" )
1166  {
1167  cout<<"zerophase"<<endl;
1168 
1169  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1170  + filePost +" "+ filestreak +" "+ fileinterp + " "+ alpchar;
1171 
1172  cout<<syscall.c_str()<<endl;
1173  if(system(syscall.c_str()))
1174  {
1175  ASSERTL0(false,syscall.c_str());
1176  }
1177 
1178  //move the advPost mesh (remark update alpha!!!)
1179  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1180  + filePost + " " + filestreak + " " + filePost + " "+ alpchar;
1181  cout<<syscall.c_str()<<endl;
1182  if(system(syscall.c_str()))
1183  {
1184  ASSERTL0(false,syscall.c_str());
1185  }
1186 
1187 
1188 
1189  //save oldstreak
1190  string oldstreak = m_sessionName +"_streak_"+ c +".fld";
1191  syscall = "cp -f " + filestreak + " " + oldstreak;
1192  cout<<syscall.c_str()<<endl;
1193  if(system(syscall.c_str()))
1194  {
1195  ASSERTL0(false,syscall.c_str());
1196  }
1197 
1198  //interpolate the streak field into the new mesh
1199  string movedmesh = m_sessionName + "_advPost_moved.xml";
1200  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1201 
1202  //create the interp streak
1203 
1204  syscall = "../../utilities/PostProcessing/Extras/FieldToField "
1205  + fileinterp + " " + filestreak + " " + movedinterpmesh + " "
1206  + interpstreak;
1207 
1208  cout<<syscall.c_str()<<endl;
1209  if(system(syscall.c_str()))
1210  {
1211  ASSERTL0(false,syscall.c_str());
1212  }
1213 
1214 
1215  //save the old mesh
1216  string meshfile = m_sessionName + ".xml";
1217  string meshold = m_sessionName +"_"+ c +".xml";
1218  syscall = "cp -f " + meshfile + " " + meshold;
1219  cout<<syscall.c_str()<<endl;
1220  if(system(syscall.c_str()))
1221  {
1222  ASSERTL0(false,syscall.c_str());
1223  }
1224 
1225  //overwriting the meshfile with the new mesh
1226  syscall = "cp -f " + movedmesh + " " + meshfile;
1227  cout<<syscall.c_str()<<endl;
1228  if(system(syscall.c_str()))
1229  {
1230  ASSERTL0(false,syscall.c_str());
1231  }
1232 
1233  //overwriting the streak file!!
1234  syscall = "cp -f " + interpstreak + " " + filestreak;
1235  cout<<syscall.c_str()<<endl;
1236  if(system(syscall.c_str()))
1237  {
1238  ASSERTL0(false,syscall.c_str());
1239  }
1240 
1241  //calculate the wave
1242  ExecuteWave();
1243 
1244  //save the wave field:
1245  string oldwave = m_sessionName +"_wave_"+c +".fld";
1246  syscall = "cp -f " + m_sessionName+".fld" + " " + oldwave;
1247  cout<<syscall.c_str()<<endl;
1248  if(system(syscall.c_str()))
1249  {
1250  ASSERTL0(false,syscall.c_str());
1251  }
1252 
1253  //save old jump conditions:
1254  string ujump = m_sessionName+"_u_5.bc";
1255  syscall = "cp -f " + ujump + " " + m_sessionName+"_u_5.bc_"+c;
1256  cout<<syscall.c_str()<<endl;
1257  if(system(syscall.c_str()))
1258  {
1259  ASSERTL0(false,syscall.c_str());
1260  }
1261 
1262  string vjump = m_sessionName+"_v_5.bc";
1263  syscall = "cp -f " + vjump + " " + m_sessionName+"_v_5.bc_"+c;
1264  cout<<syscall.c_str()<<endl;
1265  if(system(syscall.c_str()))
1266  {
1267  ASSERTL0(false,syscall.c_str());
1268  }
1269  cnt++;
1270 
1271 
1272 
1273 
1274  //use relaxation
1275  if(GetVWIIterationType()!=eFixedWaveForcingWithSubIterationOnAlpha
1276  )
1277  {
1278  // the critical layer should be the bnd region 3
1279  //int reg =3;
1280  //FileRelaxation(reg);
1281  }
1282  char c1[16]="";
1283  sprintf(c1,"%d",cnt);
1284  //calculate the jump conditions
1285  string wavefile = m_sessionName +".fld";
1286  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs "
1287  + movedmesh + " " + wavefile + " " + interpstreak + "> data"+c1;
1288  cout<<syscall.c_str()<<endl;
1289  if(system(syscall.c_str()))
1290  {
1291  ASSERTL0(false,syscall.c_str());
1292  }
1293 
1294  //move the new name_interp_moved.xml into name_interp.xml
1295  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1296  cout<<syscall.c_str()<<endl;
1297  if(system(syscall.c_str()))
1298  {
1299  ASSERTL0(false,syscall.c_str());
1300  }
1301  //move the new name_advPost_moved.xml into name_advPost.xml
1302  syscall = "cp -f " + movedmesh + " " + filePost;
1303  cout<<syscall.c_str()<<endl;
1304  if(system(syscall.c_str()))
1305  {
1306  ASSERTL0(false,syscall.c_str());
1307  }
1308 
1309 
1310  }
1311  else if( m_sessionVWI->GetSolverInfo("INTERFACE")=="phase" )
1312  {
1313 cout<<"phase"<<endl;
1314  //determine cr:
1315  NekDouble cr;
1316  string cr_str;
1317  stringstream st;
1318 
1319  //calculate the wave
1320  ExecuteWave();
1321 
1322  //save oldstreak
1323  string oldstreak = m_sessionName +"_streak_"+ c +".fld";
1324  syscall = "cp -f " + filestreak + " " + oldstreak;
1325  cout<<syscall.c_str()<<endl;
1326  if(system(syscall.c_str()))
1327  {
1328  ASSERTL0(false,syscall.c_str());
1329  }
1330 
1331  //save wave
1332  syscall = "cp -f " + m_sessionName+".fld" + " " + filewave;
1333  cout<<syscall.c_str()<<endl;
1334  if(system(syscall.c_str()))
1335  {
1336  ASSERTL0(false,syscall.c_str());
1337  }
1338  //save the old mesh
1339  string meshfile = m_sessionName + ".xml";
1340  string meshold = m_sessionName +"_"+ c +".xml";
1341  syscall = "cp -f " + meshfile + " " + meshold;
1342  cout<<syscall.c_str()<<endl;
1343  if(system(syscall.c_str()))
1344  {
1345  ASSERTL0(false,syscall.c_str());
1346  }
1347 
1348  //save the oldwave field:
1349  string oldwave = m_sessionName +"_wave_"+c +".fld";
1350  syscall = "cp -f " + m_sessionName+".fld" + " " + oldwave;
1351  cout<<syscall.c_str()<<endl;
1352  if(system(syscall.c_str()))
1353  {
1354  ASSERTL0(false,syscall.c_str());
1355  }
1356 
1357  //save old jump conditions:
1358  string ujump = m_sessionName+"_u_5.bc";
1359  syscall = "cp -f " + ujump + " " + m_sessionName+"_u_5.bc_"+c;
1360  cout<<syscall.c_str()<<endl;
1361  if(system(syscall.c_str()))
1362  {
1363  ASSERTL0(false,syscall.c_str());
1364  }
1365 
1366  string vjump = m_sessionName+"_v_5.bc";
1367  syscall = "cp -f " + vjump + " " + m_sessionName+"_v_5.bc_"+c;
1368  cout<<syscall.c_str()<<endl;
1369  if(system(syscall.c_str()))
1370  {
1371  ASSERTL0(false,syscall.c_str());
1372  }
1373  cnt++;
1374 
1375 
1376  cr = m_leading_imag_evl[0]/m_alpha[0];
1377  st << cr;
1378  cr_str = st.str();
1379 cout<<"cr="<<cr_str<<endl;
1380  //NB -g or NOT!!!
1381  //move the mesh around the critical layer
1382  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1383  + filePost +" "+ filestreak +" "+ fileinterp + " "+ alpchar
1384  +" "+cr_str;
1385 
1386  cout<<syscall.c_str()<<endl;
1387  if(system(syscall.c_str()))
1388  {
1389  ASSERTL0(false,syscall.c_str());
1390  }
1391  //NB -g or NOT!!!
1392  //move the advPost mesh (remark update alpha!!!)
1393  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1394  + filePost + " " + filestreak + " " + filePost + " "+ alpchar
1395  +" "+cr_str;
1396  cout<<syscall.c_str()<<endl;
1397  if(system(syscall.c_str()))
1398  {
1399  ASSERTL0(false,syscall.c_str());
1400  }
1401 
1402  //interp streak into the new mesh
1403  syscall = "../../utilities/PostProcessing/Extras/FieldToField "
1404  + fileinterp + " " + filestreak + " " + movedinterpmesh + " "
1405  + interpstreak;
1406 
1407  cout<<syscall.c_str()<<endl;
1408  if(system(syscall.c_str()))
1409  {
1410  ASSERTL0(false,syscall.c_str());
1411  }
1412 
1413  //split wave sol
1414  syscall = "../../utilities/PostProcessing/Extras/SplitFld "
1415  + filePost + " " + filewave;
1416 
1417  cout<<syscall.c_str()<<endl;
1418  if(system(syscall.c_str()))
1419  {
1420  ASSERTL0(false,syscall.c_str());
1421  }
1422  //interp wave
1423  syscall = "../../utilities/PostProcessing/Extras/FieldToField "
1424  + filePost + " " + filewavepressure + " " + movedmesh + " "
1425  + interwavepressure;
1426 
1427  cout<<syscall.c_str()<<endl;
1428  if(system(syscall.c_str()))
1429  {
1430  ASSERTL0(false,syscall.c_str());
1431  }
1432 
1433 
1434 
1435 
1436  //use relaxation
1437  if(GetVWIIterationType()!=eFixedWaveForcingWithSubIterationOnAlpha
1438  )
1439  {
1440  // the critical layer should be the bnd region 3
1441  //int reg =3;
1442  //FileRelaxation(reg);
1443  }
1444  char c1[16]="";
1445  sprintf(c1,"%d",cnt);
1446 
1447  //cp wavepressure to m_sessionName.fld(to get
1448  // the right bcs names using FldCalcBCs)
1449  syscall = "cp -f "+ interwavepressure +" "+m_sessionName+".fld";
1450  cout<<syscall.c_str()<<endl;
1451  if(system(syscall.c_str()))
1452  {
1453  ASSERTL0(false,syscall.c_str());
1454  }
1455 
1456  //calculate the jump conditions
1457  //NB -g or NOT!!!
1458  syscall = "../../utilities/PostProcessing/Extras/FldCalcBCs "
1459  + movedmesh + " " +m_sessionName+".fld" + " "
1460  + interpstreak + "> data"+c1;
1461  cout<<syscall.c_str()<<endl;
1462  if(system(syscall.c_str()))
1463  {
1464  ASSERTL0(false,syscall.c_str());
1465  }
1466 
1467 
1468  //overwriting the meshfile with the new mesh
1469  syscall = "cp -f " + movedmesh + " " + meshfile;
1470  cout<<syscall.c_str()<<endl;
1471  if(system(syscall.c_str()))
1472  {
1473  ASSERTL0(false,syscall.c_str());
1474  }
1475 
1476  //overwriting the streak file!! (maybe is useless)
1477  syscall = "cp -f " + interpstreak + " " + filestreak;
1478  cout<<syscall.c_str()<<endl;
1479  if(system(syscall.c_str()))
1480  {
1481  ASSERTL0(false,syscall.c_str());
1482  }
1483  //move the new name_interp_moved.xml into name_interp.xml
1484  syscall = "cp -f " + movedinterpmesh + " " + fileinterp;
1485  cout<<syscall.c_str()<<endl;
1486  if(system(syscall.c_str()))
1487  {
1488  ASSERTL0(false,syscall.c_str());
1489  }
1490  //move the new name_advPost_moved.xml into name_advPost.xml
1491  syscall = "cp -f " + movedmesh + " " + filePost;
1492  cout<<syscall.c_str()<<endl;
1493  if(system(syscall.c_str()))
1494  {
1495  ASSERTL0(false,syscall.c_str());
1496  }
1497  }
1498  }
1499  else
1500  {
1501  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::EnableManagement("GlobalLinSys");
1502  ExecuteRoll();
1503  LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,MultiRegions::GlobalLinSys>::DisableManagement("GlobalLinSys");
1504 
1505 #ifndef _WIN32
1506  sleep(3);
1507 #endif
1508  ExecuteStreak();
1509 #ifndef _WIN32
1510  sleep(3);
1511 #endif
1512 
1513  if(m_moveMeshToCriticalLayer)
1514  {
1515  string syscall;
1516  char alpchar[16]="";
1517  sprintf(alpchar, "%f", m_alpha[0]);
1518 
1519  string filePost = m_sessionName + "_advPost.xml";
1520  string filestreak = m_sessionName + "_streak.fld";
1521  string filewave = m_sessionName + "_wave.fld";
1522  string filewavepressure = m_sessionName + "_wave_p_split.fld";
1523  string fileinterp = m_sessionName + "_interp.xml";
1524  string interpstreak = m_sessionName +"_interpstreak.fld";
1525  string interwavepressure = m_sessionName +"_wave_p_split_interp.fld";
1526  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1527  + filePost +" "+ filestreak +" "+ fileinterp + " "+ alpchar;
1528 
1529  cout<<syscall.c_str()<<endl;
1530  if(system(syscall.c_str()))
1531  {
1532  ASSERTL0(false,syscall.c_str());
1533  }
1534 
1535  //move the advPost mesh (remark update alpha!!!)
1536  syscall = "../../utilities/PostProcessing/Extras/MoveMesh "
1537  + filePost + " " + filestreak + " " + filePost + " "+ alpchar;
1538  cout<<syscall.c_str()<<endl;
1539  if(system(syscall.c_str()))
1540  {
1541  ASSERTL0(false,syscall.c_str());
1542  }
1543 
1544  //save oldstreak
1545  string oldstreak = m_sessionName +"_streak.fld";
1546  syscall = "cp -f " + filestreak + " " + oldstreak;
1547  cout<<syscall.c_str()<<endl;
1548  if(system(syscall.c_str()))
1549  {
1550  ASSERTL0(false,syscall.c_str());
1551  }
1552 
1553  //interpolate the streak field into the new mesh
1554  string movedmesh = m_sessionName + "_advPost_moved.xml";
1555  string movedinterpmesh = m_sessionName + "_interp_moved.xml";
1556 
1557  //create the interp streak
1558 
1559  syscall = "../../utilities/PostProcessing/Extras/FieldToField "
1560  + fileinterp + " " + filestreak + " " + movedinterpmesh
1561  + " " + interpstreak;
1562 
1563  cout<<syscall.c_str()<<endl;
1564  if(system(syscall.c_str()))
1565  {
1566  ASSERTL0(false,syscall.c_str());
1567  }
1568 
1569  //save the old mesh
1570  string meshfile = m_sessionName + ".xml";
1571  string meshold = m_sessionName + ".xml";
1572  syscall = "cp -f " + meshfile + " " + meshold;
1573  cout<<syscall.c_str()<<endl;
1574  if(system(syscall.c_str()))
1575  {
1576  ASSERTL0(false,syscall.c_str());
1577  }
1578 
1579  //overwriting the meshfile with the new mesh
1580  syscall = "cp -f " + movedmesh + " " + meshfile;
1581  cout<<syscall.c_str()<<endl;
1582  if(system(syscall.c_str()))
1583  {
1584  ASSERTL0(false,syscall.c_str());
1585  }
1586 
1587  //overwriting the streak file!!
1588  syscall = "cp -f " + interpstreak + " " + filestreak;
1589  cout<<syscall.c_str()<<endl;
1590  if(system(syscall.c_str()))
1591  {
1592  ASSERTL0(false,syscall.c_str());
1593  }
1594  }
1595 
1596  ExecuteWave();
1597 
1598  if(CalcWaveForce)
1599  {
1600  CalcNonLinearWaveForce();
1601  }
1602  }
1603  }
Nektar::LibUtilities::NekManager::DisableManagement
static void DisableManagement(std::string whichPool="")
Definition: NekManager.hpp:293
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Nektar::VortexWaveInteraction::m_sessionRoll
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Definition: VortexWaveInteraction.h:272
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Nektar::VortexWaveInteraction::GetVWIIterationType
VWIIterationType GetVWIIterationType(void)
Definition: VortexWaveInteraction.h:115
Nektar::LibUtilities::NekManager::EnableManagement
static void EnableManagement(std::string whichPool="")
Definition: NekManager.hpp:274
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
void Nektar::VortexWaveInteraction::ExecuteRoll ( void  )

Definition at line 353 of file VortexWaveInteraction.cpp.

References Nektar::IncNavierStokes::AddForcing(), CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetEquationSystemFactory(), Nektar::SolverUtils::GetForcingFactory(), m_moveMeshToCriticalLayer, m_rollField, m_rollForceScale, m_sessionName, m_sessionRoll, m_solverRoll, m_vwiForcing, m_vwiForcingObj, Vmath::Smul(), and Vmath::Vcopy().

Referenced by ExecuteLoop().

354  {
355  //set up the equation system to update the mesh
356  if(m_sessionRoll->DefinesSolverInfo("INTERFACE"))
357  {
358  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
359  EquationSystemSharedPtr solverRoll = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionRoll);
360  // The forcing terms are inserted as N bcs
361  // Execute Roll
362  cout << "Executing Roll solver" << endl;
363  solverRoll->DoInitialise();
364  solverRoll->DoSolve();
365  solverRoll->Output();
366  m_rollField = solverRoll->UpdateFields();
367  for(int g=0; g< solverRoll->GetNvariables(); ++g)
368  {
369  NekDouble vL2Error = solverRoll->L2Error(g,false);
370  NekDouble vLinfError = solverRoll->LinfError(g);
371  cout << "L 2 error (variable " << solverRoll->GetVariable(g) << ") : " << vL2Error << endl;
372  cout << "L inf error (variable " << solverRoll->GetVariable(g) << ") : " << vLinfError << endl;
373  }
374  }
375  else
376  {
377  if(m_moveMeshToCriticalLayer)
378  {
379  string vEquation = m_sessionRoll->GetSolverInfo("solvertype");
380  EquationSystemSharedPtr solverRoll = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionRoll);
381  }
382  else
383  {
384  const int npoints = m_solverRoll->GetNpoints();
385  const int ncoeffs = m_solverRoll->GetNcoeffs();
386 
387  static int init = 1;
388  if(init)
389  {
390  m_solverRoll->DoInitialise();
391  m_vwiForcingObj = boost::dynamic_pointer_cast<SolverUtils::ForcingProgrammatic>(GetForcingFactory().CreateInstance("Programmatic", m_sessionRoll, m_solverRoll->UpdateFields(), m_solverRoll->UpdateFields().num_elements() - 1, 0));
392 
393  std::vector<std::string> vFieldNames = m_sessionRoll->GetVariables();
394  vFieldNames.erase(vFieldNames.end()-1);
395 
396  m_solverRoll->EvaluateFunction(vFieldNames, m_vwiForcingObj->UpdateForces(), "BodyForce");
397 
398  // Scale forcing
399  for(int i = 0; i < m_vwiForcingObj->UpdateForces().num_elements(); ++i)
400  {
401  m_solverRoll->UpdateFields()[0]->FwdTrans(m_vwiForcingObj->UpdateForces()[i], m_vwiForcing[2+i]);
402  Vmath::Smul(npoints,m_rollForceScale,m_vwiForcingObj->UpdateForces()[i],1,m_vwiForcingObj->UpdateForces()[i],1);
403  }
404 
405  IncNavierStokesSharedPtr ins = m_solverRoll->as<IncNavierStokes>();
406  ins->AddForcing(m_vwiForcingObj);
407 
408  init = 0;
409  }
410  else // use internal definition of forcing in m_vwiForcing
411  {
412  // Scale forcing
413  for(int i = 0; i < m_vwiForcingObj->UpdateForces().num_elements(); ++i)
414  {
415  m_solverRoll->UpdateFields()[i]->BwdTrans(m_vwiForcing[i],m_vwiForcingObj->UpdateForces()[i]);
416  Vmath::Smul(npoints,m_rollForceScale,m_vwiForcingObj->UpdateForces()[i],1,m_vwiForcingObj->UpdateForces()[i],1);
417  }
418 
419  // Shift m_vwiForcing for new restart in case of relaxation
420  Vmath::Vcopy(ncoeffs,m_vwiForcing[0],1,m_vwiForcing[2],1);
421  Vmath::Vcopy(ncoeffs,m_vwiForcing[1],1,m_vwiForcing[3],1);
422  }
423  }
424 
425  // Execute Roll
426  cout << "Executing Roll solver" << endl;
427  m_solverRoll->DoSolve();
428  m_solverRoll->Output();
429  m_rollField = m_solverRoll->UpdateFields();
430  for(int g=0; g< m_solverRoll->GetNvariables(); ++g)
431  {
432  NekDouble vL2Error = m_solverRoll->L2Error(g,false);
433  NekDouble vLinfError = m_solverRoll->LinfError(g);
434  cout << "L 2 error (variable " << m_solverRoll->GetVariable(g) << ") : " << vL2Error << endl;
435  cout << "L inf error (variable " << m_solverRoll->GetVariable(g) << ") : " << vLinfError << endl;
436  }
437 
438 
439  }
440 
441  // Copy .fld file to .rst and base.fld
442  cout << "Executing cp -f session.fld session.rst" << endl;
443  CopyFile(".fld",".rst");
444 
445  // Write out data into base flow with variable Vx,Vy
446  cout << "Writing data to session-Base.fld" << endl;
447 
448  std::vector<std::string> variables(2);
449  variables[0] = "Vx"; variables[1] = "Vy";
450  std::vector<Array<OneD, NekDouble> > outfield(2);
451  outfield[0] = m_solverRoll->UpdateFields()[0]->UpdateCoeffs();
452  outfield[1] = m_solverRoll->UpdateFields()[1]->UpdateCoeffs();
453  std::string outname = m_sessionName + "-Base.fld";
454  m_solverRoll->WriteFld(outname, m_solverRoll->UpdateFields()[0],
455  outfield, variables);
456  }
Nektar::VortexWaveInteraction::m_sessionRoll
LibUtilities::SessionReaderSharedPtr m_sessionRoll
Definition: VortexWaveInteraction.h:272
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Nektar::IncNavierStokesSharedPtr
boost::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
Definition: IncNavierStokes.h:238
Nektar::VortexWaveInteraction::m_rollField
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
Definition: VortexWaveInteraction.h:267
Nektar::SolverUtils::GetForcingFactory
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
Definition: Forcing.cpp:44
Nektar::VortexWaveInteraction::m_solverRoll
EquationSystemSharedPtr m_solverRoll
Definition: VortexWaveInteraction.h:273
Nektar::VortexWaveInteraction::m_vwiForcingObj
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
Definition: VortexWaveInteraction.h:261
Nektar::SolverUtils::EquationSystemSharedPtr
boost::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:56
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:83
Nektar::VortexWaveInteraction::m_vwiForcing
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
Definition: VortexWaveInteraction.h:260
Nektar::VortexWaveInteraction::CopyFile
void CopyFile(std::string file1end, std::string file2end)
Definition: VortexWaveInteraction.cpp:1050
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::VortexWaveInteraction::ExecuteStreak ( void  )

Definition at line 459 of file VortexWaveInteraction.cpp.

References CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetDriverFactory(), Nektar::SolverUtils::GetEquationSystemFactory(), m_sessionStreak, m_sessionVWI, and m_streakField.

Referenced by ExecuteLoop().

460  {
461  // Create driver
462 #if 1
463  std::string vDriverModule;
464  m_sessionStreak->LoadSolverInfo("Driver", vDriverModule, "Standard");
465 
466  DriverSharedPtr solverStreak = GetDriverFactory().CreateInstance(vDriverModule, m_sessionStreak);
467  solverStreak->Execute();
468 
469  m_streakField = solverStreak->GetEqu()[0]->UpdateFields();
470 #else
471  // Setup and execute Advection Diffusion solver
472  string vEquation = m_sessionStreak->GetSolverInfo("EqType");
473  EquationSystemSharedPtr solverStreak = GetEquationSystemFactory().CreateInstance(vEquation,m_sessionStreak);
474 
475  cout << "Executing Streak Solver" << endl;
476  solverStreak->DoInitialise();
477  solverStreak->DoSolve();
478  solverStreak->Output();
479 
480  m_streakField = solverStreak->UpdateFields();
481 
482  if(m_sessionVWI->DefinesSolverInfo("INTERFACE"))
483  {
484  for(int g=0; g< solverStreak->GetNvariables(); ++g)
485  {
486  NekDouble vL2Error = solverStreak->L2Error(g,false);
487  NekDouble vLinfError = solverStreak->LinfError(g);
488  cout << "L 2 error (variable " << solverStreak->GetVariable(g) << ") : " << vL2Error << endl;
489  cout << "L inf error (variable " << solverStreak->GetVariable(g) << ") : " << vLinfError << endl;
490  }
491  }
492 #endif
493 
494  cout << "Executing cp -f session.fld session_streak.fld" << endl;
495  CopyFile(".fld","_streak.fld");
496  }
Nektar::SolverUtils::DriverSharedPtr
boost::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:52
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Nektar::SolverUtils::EquationSystemSharedPtr
boost::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Definition: EquationSystem.h:56
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:83
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
Nektar::VortexWaveInteraction::CopyFile
void CopyFile(std::string file1end, std::string file2end)
Definition: VortexWaveInteraction.cpp:1050
Nektar::SolverUtils::GetDriverFactory
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:66
Nektar::VortexWaveInteraction::m_streakField
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Definition: VortexWaveInteraction.h:265
Nektar::VortexWaveInteraction::m_sessionStreak
LibUtilities::SessionReaderSharedPtr m_sessionStreak
Definition: VortexWaveInteraction.h:275
void Nektar::VortexWaveInteraction::ExecuteWave ( void  )

Do linearised NavierStokes Session with Modified Arnoldi

Definition at line 498 of file VortexWaveInteraction.cpp.

References CopyFile(), Nektar::LibUtilities::NekFactory< tKey, tBase, >::CreateInstance(), Nektar::SolverUtils::GetDriverFactory(), m_alpha, m_leading_imag_evl, m_leading_real_evl, m_sessionVWI, m_sessionWave, m_wavePressure, and m_waveVelocities.

Referenced by CalcL2ToLinfPressure(), DoFixedForcingIteration(), ExecuteLoop(), and main().

499  {
500 
501  // Set the initial beta value in stability to be equal to VWI file
502  std::string LZstr("LZ");
503  NekDouble LZ = 2*M_PI/m_alpha[0];
504  cout << "Setting LZ in Linearised solver to " << LZ << endl;
505  m_sessionWave->SetParameter(LZstr,LZ);
506 
507  // Create driver
508  std::string vDriverModule;
509  m_sessionWave->LoadSolverInfo("Driver", vDriverModule, "ModifiedArnoldi");
510  cout << "Setting up linearised NS sovler" << endl;
511  DriverSharedPtr solverWave = GetDriverFactory().CreateInstance(vDriverModule, m_sessionWave);
512 
513  /// Do linearised NavierStokes Session with Modified Arnoldi
514  cout << "Executing wave solution " << endl;
515  solverWave->Execute();
516 
517  // Copy file to a rst location for next restart
518  cout << "Executing cp -f session_eig_0 session_eig_0.rst" << endl;
519  CopyFile("_eig_0","_eig_0.rst");
520 
521  // Store data relevant to other operations
522  m_leading_real_evl[0] = solverWave->GetRealEvl()[0];
523  m_leading_imag_evl[0] = solverWave->GetImagEvl()[0];
524 
525  // note this will only be true for modified Arnoldi
526  NekDouble realShift = 0.0;
527  if(m_sessionWave->DefinesParameter("RealShift"))
528  {
529  bool defineshift;
530  // only use shift in modifiedArnoldi solver since
531  // implicitly handled in Arpack.
532  m_sessionWave->MatchSolverInfo("Driver","ModifiedArnoldi",defineshift,true);
533  if(defineshift)
534  {
535  realShift = m_sessionWave->GetParameter("RealShift");
536  }
537 
538  }
539 
540  // Set up leading eigenvalue from inverse
541  NekDouble invmag = 1.0/(m_leading_real_evl[0]*m_leading_real_evl[0] +
542  m_leading_imag_evl[0]*m_leading_imag_evl[0]);
543  m_leading_real_evl[0] *= -invmag;
544  m_leading_real_evl[0] += realShift;
545  m_leading_imag_evl[0] *= invmag;
546 
547 
548  m_waveVelocities = solverWave->GetEqu()[0]->UpdateFields();
549  m_wavePressure = solverWave->GetEqu()[0]->GetPressure();
550 
551 
552  if( m_sessionVWI->DefinesSolverInfo("INTERFACE") )
553  {
554  cout << "Growth =" <<m_leading_real_evl[0]<<endl;
555  cout << "Phase =" <<m_leading_imag_evl[0]<<endl;
556  }
557 
558  }
Nektar::VortexWaveInteraction::m_wavePressure
MultiRegions::ExpListSharedPtr m_wavePressure
Definition: VortexWaveInteraction.h:258
Nektar::SolverUtils::DriverSharedPtr
boost::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:52
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Nektar::VortexWaveInteraction::m_waveVelocities
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
Definition: VortexWaveInteraction.h:257
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:241
Nektar::VortexWaveInteraction::m_sessionWave
LibUtilities::SessionReaderSharedPtr m_sessionWave
Definition: VortexWaveInteraction.h:276
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
Nektar::VortexWaveInteraction::CopyFile
void CopyFile(std::string file1end, std::string file2end)
Definition: VortexWaveInteraction.cpp:1050
Nektar::SolverUtils::GetDriverFactory
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:66
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
void Nektar::VortexWaveInteraction::FileRelaxation ( int  reg)

====================================================

Todo:

Please update to use MeshGraph::Read(vSession)

Definition at line 2005 of file VortexWaveInteraction.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), m_bcsForcing, m_rollField, m_sessionName, m_sessionVWI, m_vwiRelaxation, Nektar::SpatialDomains::MeshGraph::Read(), Vmath::Smul(), Vmath::Svtvp(), and Vmath::Vcopy().

2006  {
2007  cout<<"relaxation..."<<endl;
2008  static int cnt=0;
2009  Array<OneD, MultiRegions::ExpListSharedPtr> Iexp
2010  =m_rollField[0]->GetBndCondExpansions();
2011  //cast to 1D explist (otherwise appenddata doesn't work)
2012  MultiRegions::ExpList1DSharedPtr Ilayer;
2013  Ilayer = MemoryManager<MultiRegions::ExpList1D>::
2014  AllocateSharedPtr(
2015  *boost::static_pointer_cast<MultiRegions::ExpList1D>(Iexp[reg]));
2016  int nq = Ilayer->GetTotPoints();
2017  if( cnt==0)
2018  {
2019  m_bcsForcing = Array<OneD, Array<OneD, NekDouble> > (4);
2020  m_bcsForcing[0] = Array<OneD, NekDouble> (4*nq);
2021  for(int i = 1; i < 4; ++i)
2022  {
2023  m_bcsForcing[i] = m_bcsForcing[i-1] + nq;
2024  }
2025  }
2026 
2027  // Read in mesh from input file
2028 
2029  /// ====================================================
2030  /// \todo Please update to use MeshGraph::Read(vSession)
2031  /// ====================================================
2032  SpatialDomains::MeshGraphSharedPtr graphShPt =
2033  SpatialDomains::MeshGraph::Read(m_sessionName+".xml");
2034 
2035 
2036  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
2037  std::vector<std::vector<NekDouble> > FieldData_u;
2038  string file = m_sessionName;
2039 
2040 
2041  file += "_u_5.bc";
2042  LibUtilities::FieldIOSharedPtr fld =
2043  MemoryManager<LibUtilities::FieldIO>::AllocateSharedPtr(m_sessionVWI->GetComm());
2044  fld->Import(file,FieldDef_u, FieldData_u);
2045  Ilayer->ExtractDataToCoeffs(FieldDef_u[0], FieldData_u[0], FieldDef_u[0]->m_fields[0],Ilayer->UpdateCoeffs());
2046  Ilayer->BwdTrans_IterPerExp(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2047 
2048  if(cnt==0)
2049  {
2050  Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[2],1);
2051  }
2052  Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[0],1);
2053 
2054  //case relaxation==0 an additional array is needed
2055  Array<OneD, NekDouble> tmp_forcing(nq, 0.0);
2056 
2057  if(cnt!=0)
2058  {
2059  if(m_vwiRelaxation==1.0)
2060  {
2061  Vmath::Vcopy(nq, m_bcsForcing[0],1, tmp_forcing,1);
2062  }
2063  Vmath::Smul(nq,1.0-m_vwiRelaxation,
2064  m_bcsForcing[0],1,m_bcsForcing[0],1);
2065  Vmath::Svtvp(nq,m_vwiRelaxation,m_bcsForcing[2],1,
2066  m_bcsForcing[0],1,Ilayer->UpdatePhys(),1);
2067  //generate again the bcs files:
2068 
2069  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);
2070  Ilayer->FwdTrans_IterPerExp(Ilayer->GetPhys(),Ilayer->UpdateCoeffs());
2071  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2072  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef1 = Ilayer->GetFieldDefinitions();
2073  std::vector<std::vector<NekDouble> > FieldData_1(FieldDef1.size());;
2074  FieldDef1[0]->m_fields.push_back("u");
2075  Ilayer->AppendFieldData(FieldDef1[0], FieldData_1[0]);
2076  fld->Write(file,FieldDef1,FieldData_1);
2077  //save the bcs for the next iteration
2078  if(m_vwiRelaxation!=1.0)
2079  {
2080  Vmath::Smul(nq,1./(1.0-m_vwiRelaxation),
2081  m_bcsForcing[0],1,m_bcsForcing[0],1);
2082  Vmath::Vcopy(nq,m_bcsForcing[0],1,m_bcsForcing[2],1);
2083  }
2084  else
2085  {
2086  Vmath::Vcopy(nq, tmp_forcing,1, m_bcsForcing[2],1);
2087  }
2088  }
2089 
2090 
2091 
2092  file = m_sessionName+ "_v_5.bc";
2093 
2094  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_v;
2095  std::vector<std::vector<NekDouble> > FieldData_v;
2096  fld->Import(file,FieldDef_v, FieldData_v);
2097  Ilayer->ExtractDataToCoeffs(FieldDef_v[0], FieldData_v[0], FieldDef_v[0]->m_fields[0],Ilayer->UpdateCoeffs());
2098  Ilayer->BwdTrans_IterPerExp(Ilayer->GetCoeffs(), Ilayer->UpdatePhys());
2099  if(cnt==0)
2100  {
2101  Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[3],1);
2102  }
2103  Vmath::Vcopy(nq,Ilayer->UpdatePhys(),1,m_bcsForcing[1],1);
2104  if(cnt!=0)
2105  {
2106  if(m_vwiRelaxation==1.0)
2107  {
2108  Vmath::Vcopy(nq, m_bcsForcing[1],1, tmp_forcing,1);
2109  }
2110  Vmath::Smul(nq,1.0-m_vwiRelaxation,
2111  m_bcsForcing[1],1,m_bcsForcing[1],1);
2112  Vmath::Svtvp(nq,m_vwiRelaxation,m_bcsForcing[3],1,
2113  m_bcsForcing[1],1,Ilayer->UpdatePhys(),1);
2114  //generate again the bcs files:
2115  Array<OneD, Array<OneD, NekDouble> > fieldcoeffs(1);
2116  Ilayer->FwdTrans_IterPerExp(Ilayer->GetPhys(),Ilayer->UpdateCoeffs());
2117  fieldcoeffs[0] = Ilayer->UpdateCoeffs();
2118  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef2 = Ilayer->GetFieldDefinitions();
2119  std::vector<std::vector<NekDouble> > FieldData_2(FieldDef2.size());;
2120  FieldDef2[0]->m_fields.push_back("v");
2121  Ilayer->AppendFieldData(FieldDef2[0], FieldData_2[0]);
2122  fld->Write(file,FieldDef2,FieldData_2);
2123  //save the bcs for the next iteration
2124  if(m_vwiRelaxation!=1.0)
2125  {
2126  Vmath::Smul(nq,1./(1.0-m_vwiRelaxation),
2127  m_bcsForcing[1],1,m_bcsForcing[1],1);
2128  Vmath::Vcopy(nq,m_bcsForcing[1],1,m_bcsForcing[3],1);
2129  }
2130  else
2131  {
2132  Vmath::Vcopy(nq, tmp_forcing,1, m_bcsForcing[3],1);
2133  }
2134 
2135 
2136  }
2137 
2138 
2139  cnt++;
2140  }
Nektar::SpatialDomains::MeshGraph::Read
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:121
Nektar::MemoryManager::AllocateSharedPtr
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:235
Nektar::VortexWaveInteraction::m_rollField
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
Definition: VortexWaveInteraction.h:267
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
Nektar::VortexWaveInteraction::m_bcsForcing
Array< OneD, Array< OneD, NekDouble > > m_bcsForcing
Definition: VortexWaveInteraction.h:263
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::MultiRegions::ExpList1DSharedPtr
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
Nektar::LibUtilities::FieldIOSharedPtr
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:225
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
Nektar::SpatialDomains::MeshGraphSharedPtr
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
NekDouble Nektar::VortexWaveInteraction::GetAlpha ( void  )
inline

Definition at line 130 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration(), and main().

131  {
132  return m_alpha[0];
133  }
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
NekDouble Nektar::VortexWaveInteraction::GetAlphaStep ( void  )
inline

Definition at line 136 of file VortexWaveInteraction.h.

Referenced by main().

137  {
138  return m_alphaStep;
139  }
NekDouble Nektar::VortexWaveInteraction::GetDAlphaDWaveForceMag ( void  )
inline

Definition at line 151 of file VortexWaveInteraction.h.

152  {
153  return m_dAlphaDWaveForceMag;
154  }
NekDouble Nektar::VortexWaveInteraction::GetEigRelTol ( void  )
inline

Definition at line 162 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

164  {
165  return m_eigRelTol;
166  }
int Nektar::VortexWaveInteraction::GetIterEnd ( )
inline

Definition at line 110 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

111  {
112  return m_iterEnd;
113  }
int Nektar::VortexWaveInteraction::GetIterStart ( )
inline

Definition at line 105 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

106  {
107  return m_iterStart;
108  }
int Nektar::VortexWaveInteraction::GetMaxOuterIterations ( void  )
inline

Definition at line 125 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

126  {
127  return m_maxOuterIterations;
128  }
int Nektar::VortexWaveInteraction::GetMaxWaveForceMagIter ( void  )
inline

Definition at line 156 of file VortexWaveInteraction.h.

Referenced by main().

157  {
158  return m_maxWaveForceMagIter;
159  }
int Nektar::VortexWaveInteraction::GetMinInnerIterations ( void  )
inline

Definition at line 168 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

169  {
170  return m_minInnerIterations;
171  }
int Nektar::VortexWaveInteraction::GetNOuterIterations ( void  )
inline

Definition at line 120 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

121  {
122  return m_nOuterIterations;
123  }
NekDouble Nektar::VortexWaveInteraction::GetPrevAlpha ( void  )
inline

Definition at line 173 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

174  {
175  return m_prevAlpha;
176  }
Array< OneD, int > Nektar::VortexWaveInteraction::GetReflectionIndex ( void  )

Definition at line 1925 of file VortexWaveInteraction.cpp.

References ASSERTL0, Nektar::LibUtilities::eTriangle, m_sessionVWI, m_waveVelocities, npts, and Vmath::Vmax().

Referenced by CalcNonLinearWaveForce().

1926  {
1927  int i,j;
1928  int npts = m_waveVelocities[0]->GetPlane(0)->GetNpoints();
1929  int nel = m_waveVelocities[0]->GetNumElmts();
1930  Array<OneD, int> index(npts);
1931 
1932  Array<OneD, NekDouble> coord(2);
1933  Array<OneD, NekDouble> coord_x(npts);
1934  Array<OneD, NekDouble> coord_y(npts);
1935 
1936  //-> Dermine the point which is on coordinate (x -> -x + Lx/2, y-> -y)
1937  m_waveVelocities[0]->GetPlane(0)->GetCoords(coord_x,coord_y);
1938  NekDouble xmax = Vmath::Vmax(npts,coord_x,1);
1939  //NekDouble tol = NekConstants::kGeomFactorsTol*NekConstants::kGeomFactorsTol;
1940  NekDouble tol = 1e-5;
1941  NekDouble xnew,ynew;
1942 
1943  int start = npts-1;
1944  int e_npts;
1945 
1946  bool useOnlyQuads = false;
1947  if(m_sessionVWI->DefinesSolverInfo("SymmetriseOnlyQuads"))
1948  {
1949  useOnlyQuads = true;
1950  }
1951 
1952  int cnt;
1953  for(int e = 0; e < nel; ++e)
1954  {
1955  e_npts = m_waveVelocities[0]->GetExp(e)->GetTotPoints();
1956  cnt = m_waveVelocities[0]->GetPhys_Offset(e);
1957 
1958  if(useOnlyQuads)
1959  {
1960  if(m_waveVelocities[0]->GetExp(e)->DetShapeType() == LibUtilities::eTriangle)
1961  {
1962  for(i = 0; i < e_npts; ++i)
1963  {
1964  index[cnt+i] = -1;
1965  }
1966  continue;
1967  }
1968  }
1969 
1970  for(i = cnt; i < cnt+e_npts; ++i)
1971  {
1972  xnew = - coord_x[i] + xmax;
1973  ynew = - coord_y[i];
1974 
1975  for(j = start; j >=0 ; --j)
1976  {
1977  if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1978  {
1979  index[i] = j;
1980  start = j;
1981  break;
1982  }
1983  }
1984 
1985  if(j == -1)
1986  {
1987 
1988  for(j = npts-1; j > start; --j)
1989  {
1990 
1991  if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1992  {
1993  index[i] = j;
1994  break;
1995  }
1996  }
1997  ASSERTL0(j != start,"Failed to find matching point");
1998  }
1999  }
2000  }
2001  return index;
2002  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:765
Nektar::VortexWaveInteraction::m_waveVelocities
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
Definition: VortexWaveInteraction.h:257
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
VWIIterationType Nektar::VortexWaveInteraction::GetVWIIterationType ( void  )
inline

Definition at line 115 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration(), ExecuteLoop(), and main().

116  {
117  return m_VWIIterationType;
118  }
NekDouble Nektar::VortexWaveInteraction::GetWaveForceMag ( void  )
inline

Definition at line 141 of file VortexWaveInteraction.h.

Referenced by main().

142  {
143  return m_waveForceMag[0];
144  }
Nektar::VortexWaveInteraction::m_waveForceMag
Array< OneD, NekDouble > m_waveForceMag
Definition: VortexWaveInteraction.h:236
NekDouble Nektar::VortexWaveInteraction::GetWaveForceMagStep ( void  )
inline

Definition at line 146 of file VortexWaveInteraction.h.

Referenced by main().

147  {
148  return m_waveForceMagStep;
149  }
bool Nektar::VortexWaveInteraction::IfIterInterface ( void  )
inline

Definition at line 209 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

210  {
211  return m_iterinterface;
212  }
void Nektar::VortexWaveInteraction::MoveFile ( std::string  fileend,
std::string  dir,
int  n 
)

Definition at line 1029 of file VortexWaveInteraction.cpp.

References ASSERTL0.

Referenced by DoFixedForcingIteration().

1030  {
1031  static map<string,int> opendir;
1032 
1033  if(opendir.find(dir) == opendir.end())
1034  {
1035  // make directory and presume will fail if it already exists
1036  string mkdir = "mkdir " + dir;
1037  system(mkdir.c_str());
1038  opendir[dir] = 1;
1039  }
1040 
1041  string savefile = dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1042  string syscall = "mv -f " + file + " " + savefile;
1043 
1044  if(system(syscall.c_str()))
1045  {
1046  ASSERTL0(false,syscall.c_str());
1047  }
1048  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
void Nektar::VortexWaveInteraction::SaveFile ( std::string  fileend,
std::string  dir,
int  n 
)

Definition at line 1005 of file VortexWaveInteraction.cpp.

References ASSERTL0.

Referenced by SaveLoopDetails().

1006  {
1007  static map<string,int> opendir;
1008 
1009  if(opendir.find(dir) == opendir.end())
1010  {
1011  // make directory and presume will fail if it already exists
1012  string mkdir = "mkdir " + dir;
1013  system(mkdir.c_str());
1014 
1015  opendir[dir] = 1;
1016  }
1017 
1018  string savefile = dir + "/" + file + "." + boost::lexical_cast<std::string>(n);
1019  string syscall = "cp -f " + file + " " + savefile;
1020 
1021  if(system(syscall.c_str()))
1022  {
1023  ASSERTL0(false,syscall.c_str());
1024  }
1025 
1026  }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
void Nektar::VortexWaveInteraction::SaveLoopDetails ( std::string  dir,
int  i 
)

Definition at line 1078 of file VortexWaveInteraction.cpp.

References m_sessionName, m_sessionVWI, and SaveFile().

Referenced by DoFixedForcingIteration().

1080  {
1081  // Save NS restart file
1082  SaveFile(m_sessionName + ".rst",SaveDir,i+1);
1083  // Save Streak Solution
1084  SaveFile(m_sessionName + "_streak.fld",SaveDir,i);
1085  // Save Wave solution output
1086  SaveFile(m_sessionName + ".evl",SaveDir,i);
1087  SaveFile(m_sessionName + "_eig_0",SaveDir,i);
1088  // Save new field file of eigenvalue
1089  SaveFile(m_sessionName + ".fld",SaveDir,i);
1090  if(!(m_sessionVWI->DefinesSolverInfo("INTERFACE")))
1091  {
1092  // Save new forcing file
1093  SaveFile(m_sessionName + ".vwi",SaveDir,i+1);
1094  }
1095 
1096  }
Nektar::VortexWaveInteraction::SaveFile
void SaveFile(std::string fileend, std::string dir, int n)
Definition: VortexWaveInteraction.cpp:1005
Nektar::VortexWaveInteraction::m_sessionVWI
LibUtilities::SessionReaderSharedPtr m_sessionVWI
Definition: VortexWaveInteraction.h:270
void Nektar::VortexWaveInteraction::SetAlpha ( NekDouble  alpha)
inline

Definition at line 178 of file VortexWaveInteraction.h.

Referenced by main().

179  {
180  m_alpha[0] = alpha;
181  }
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
void Nektar::VortexWaveInteraction::SetAlphaStep ( NekDouble  step)
inline

Definition at line 194 of file VortexWaveInteraction.h.

195  {
196  m_alphaStep = step;
197  }
void Nektar::VortexWaveInteraction::SetEigRelTol ( NekDouble  tol)
inline

Definition at line 189 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

190  {
191  m_eigRelTol = tol;
192  }
void Nektar::VortexWaveInteraction::SetMinInnerIterations ( int  niter)
inline

Definition at line 199 of file VortexWaveInteraction.h.

Referenced by DoFixedForcingIteration().

200  {
201  m_minInnerIterations = niter;
202  }
void Nektar::VortexWaveInteraction::SetPrevAlpha ( NekDouble  alpha)
inline

Definition at line 204 of file VortexWaveInteraction.h.

Referenced by main().

205  {
206  m_prevAlpha = alpha;
207  }
void Nektar::VortexWaveInteraction::SetWaveForceMag ( NekDouble  mag)
inline

Definition at line 184 of file VortexWaveInteraction.h.

Referenced by main().

185  {
186  m_waveForceMag[0] = mag;
187  }
Nektar::VortexWaveInteraction::m_waveForceMag
Array< OneD, NekDouble > m_waveForceMag
Definition: VortexWaveInteraction.h:236
void Nektar::VortexWaveInteraction::UpdateAlpha ( int  n)

Definition at line 1812 of file VortexWaveInteraction.cpp.

References Vmath::Imin(), m_alpha, m_alphaStep, m_leading_imag_evl, m_leading_real_evl, and Vmath::Vcopy().

Referenced by DoFixedForcingIteration(), and VortexWaveInteraction().

1813  {
1814  NekDouble alp_new = 0.0;
1815 
1816 
1817  if(outeriter == 1)
1818  {
1819  m_alpha[1] = m_alpha[0];
1820  if(m_leading_real_evl[0] > 0.0)
1821  {
1822  alp_new = m_alpha[0] + m_alphaStep;
1823  }
1824  else
1825  {
1826  alp_new = m_alpha[0] - m_alphaStep;
1827  }
1828  }
1829  else
1830  {
1831  int i;
1832  int nstore = (m_alpha.num_elements() < outeriter)? m_alpha.num_elements(): outeriter;
1833  Array<OneD, NekDouble> Alpha(nstore);
1834  Array<OneD, NekDouble> Growth(nstore);
1835 
1836  Vmath::Vcopy(nstore,m_alpha,1,Alpha,1);
1837  Vmath::Vcopy(nstore,m_leading_real_evl,1,Growth,1);
1838 
1839  // Sort Alpha Growth values;
1840  NekDouble store;
1841  int k;
1842  for(i = 0; i < nstore; ++i)
1843  {
1844  k = Vmath::Imin(nstore-i,&Alpha[i],1);
1845 
1846  store = Alpha[i];
1847  Alpha[i] = Alpha[i+k];
1848  Alpha[i+k] = store;
1849 
1850  store = Growth[i];
1851  Growth[i] = Growth[i+k];
1852  Growth[i+k] = store;
1853  }
1854 
1855  // See if we have any values that cross zero
1856  for(i = 0; i < nstore-1; ++i)
1857  {
1858  if(Growth[i]*Growth[i+1] < 0.0)
1859  {
1860  break;
1861  }
1862  }
1863 
1864  if(i != nstore-1)
1865  {
1866  if(nstore == 2)
1867  {
1868  alp_new = (Alpha[0]*Growth[1] - Alpha[1]*Growth[0])/(Growth[1]-Growth[0]);
1869  }
1870  else
1871  {
1872  // use a quadratic fit and step through 10000 points
1873  // to find zero.
1874  int j;
1875  int nsteps = 10000;
1876  int idx = (i == 0)?1:i;
1877  NekDouble da = Alpha[idx+1] - Alpha[idx-1];
1878  NekDouble gval_m1 = Growth[idx-1],a,gval;
1879  NekDouble c1 = Growth[idx-1]/(Alpha[idx-1]-Alpha[idx])/
1880  (Alpha[idx-1]-Alpha[idx+1]);
1881  NekDouble c2 = Growth[idx]/(Alpha[idx]-Alpha[idx-1])/
1882  (Alpha[idx]-Alpha[idx+1]);
1883  NekDouble c3 = Growth[idx+1]/(Alpha[idx+1]-Alpha[idx-1])/
1884  (Alpha[idx+1]-Alpha[idx]);
1885 
1886  for(j = 1; j < nsteps+1; ++j)
1887  {
1888  a = Alpha[i] + j*da/(NekDouble) nsteps;
1889  gval = c1*(a - Alpha[idx ])*(a - Alpha[idx+1])
1890  + c2*(a - Alpha[idx-1])*(a - Alpha[idx+1])
1891  + c3*(a - Alpha[idx-1])*(a - Alpha[idx]);
1892 
1893  if(gval*gval_m1 < 0.0)
1894  {
1895  alp_new = ((a+da/(NekDouble)nsteps)*gval - a*gval_m1)/
1896  (gval - gval_m1);
1897  break;
1898  }
1899  }
1900  }
1901  }
1902  else // step backward/forward
1903  {
1904  if(Growth[i] > 0.0)
1905  {
1906  alp_new = m_alpha[0] + m_alphaStep;
1907  }
1908  else
1909  {
1910  alp_new = m_alpha[0] - m_alphaStep;
1911  }
1912  }
1913  }
1914 
1915  for(int i = m_alpha.num_elements()-1; i > 0; --i)
1916  {
1917  m_alpha[i] = m_alpha[i-1];
1918  m_leading_real_evl[i] = m_leading_real_evl[i-1];
1919  m_leading_imag_evl[i] = m_leading_imag_evl[i-1];
1920  }
1921 
1922  m_alpha[0] = alp_new;
1923  }
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:833
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:241
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
void Nektar::VortexWaveInteraction::UpdateDAlphaDWaveForceMag ( NekDouble  alphainit)

Definition at line 1807 of file VortexWaveInteraction.cpp.

References m_alpha, m_dAlphaDWaveForceMag, and m_waveForceMagStep.

Referenced by DoFixedForcingIteration().

1808  {
1809  m_dAlphaDWaveForceMag = (m_alpha[0]-alpha_init)/m_waveForceMagStep;
1810  }
Nektar::VortexWaveInteraction::m_alpha
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Definition: VortexWaveInteraction.h:244
void Nektar::VortexWaveInteraction::UpdateWaveForceMag ( int  n)

Definition at line 1693 of file VortexWaveInteraction.cpp.

References Vmath::Imin(), m_leading_imag_evl, m_leading_real_evl, m_waveForceMag, m_waveForceMagStep, and Vmath::Vcopy().

Referenced by DoFixedForcingIteration().

1694  {
1695  NekDouble wavef_new = 0.0;
1696 
1697 
1698  if(outeriter == 1)
1699  {
1700  m_waveForceMag[1] = m_waveForceMag[0];
1701  if(m_leading_real_evl[0] > 0.0)
1702  {
1703  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1704  }
1705  else
1706  {
1707  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1708  }
1709  }
1710  else
1711  {
1712  int i;
1713  int nstore = (m_waveForceMag.num_elements() < outeriter)? m_waveForceMag.num_elements(): outeriter;
1714  Array<OneD, NekDouble> WaveForce(nstore);
1715  Array<OneD, NekDouble> Growth(nstore);
1716 
1717  Vmath::Vcopy(nstore,m_waveForceMag,1,WaveForce,1);
1718  Vmath::Vcopy(nstore,m_leading_real_evl,1,Growth,1);
1719 
1720  // Sort WaveForce Growth values;
1721  NekDouble store;
1722  int k;
1723  for(i = 0; i < nstore; ++i)
1724  {
1725  k = Vmath::Imin(nstore-i,&WaveForce[i],1);
1726 
1727  store = WaveForce[i];
1728  WaveForce[i] = WaveForce[i+k];
1729  WaveForce[i+k] = store;
1730 
1731  store = Growth[i];
1732  Growth[i] = Growth[i+k];
1733  Growth[i+k] = store;
1734  }
1735 
1736  // See if we have any values that cross zero
1737  for(i = 0; i < nstore-1; ++i)
1738  {
1739  if(Growth[i]*Growth[i+1] < 0.0)
1740  {
1741  break;
1742  }
1743  }
1744 
1745  if(i != nstore-1)
1746  {
1747  if(nstore == 2)
1748  {
1749  wavef_new = (WaveForce[0]*Growth[1] - WaveForce[1]*Growth[0])/(Growth[1]-Growth[0]);
1750  }
1751  else
1752  {
1753  // use a quadratic fit and step through 10000 points
1754  // to find zero.
1755  int j;
1756  int nsteps = 10000;
1757  int idx = (i == 0)?1:i;
1758  NekDouble da = WaveForce[idx+1] - WaveForce[idx-1];
1759  NekDouble gval_m1 = Growth[idx-1],a,gval;
1760  NekDouble c1 = Growth[idx-1]/(WaveForce[idx-1]-WaveForce[idx])/
1761  (WaveForce[idx-1]-WaveForce[idx+1]);
1762  NekDouble c2 = Growth[idx]/(WaveForce[idx]-WaveForce[idx-1])/
1763  (WaveForce[idx]-WaveForce[idx+1]);
1764  NekDouble c3 = Growth[idx+1]/(WaveForce[idx+1]-WaveForce[idx-1])/
1765  (WaveForce[idx+1]-WaveForce[idx]);
1766 
1767  for(j = 1; j < nsteps+1; ++j)
1768  {
1769  a = WaveForce[i] + j*da/(NekDouble) nsteps;
1770  gval = c1*(a - WaveForce[idx ])*(a - WaveForce[idx+1])
1771  + c2*(a - WaveForce[idx-1])*(a - WaveForce[idx+1])
1772  + c3*(a - WaveForce[idx-1])*(a - WaveForce[idx]);
1773 
1774  if(gval*gval_m1 < 0.0)
1775  {
1776  wavef_new = ((a+da/(NekDouble)nsteps)*gval - a*gval_m1)/
1777  (gval - gval_m1);
1778  break;
1779  }
1780  }
1781  }
1782  }
1783  else // step backward/forward
1784  {
1785  if(Growth[i] > 0.0)
1786  {
1787  wavef_new = m_waveForceMag[0] - m_waveForceMagStep;
1788  }
1789  else
1790  {
1791  wavef_new = m_waveForceMag[0] + m_waveForceMagStep;
1792  }
1793  }
1794  }
1795 
1796  for(int i = m_waveForceMag.num_elements()-1; i > 0; --i)
1797  {
1798  m_waveForceMag[i] = m_waveForceMag[i-1];
1799  m_leading_real_evl[i] = m_leading_real_evl[i-1];
1800  m_leading_imag_evl[i] = m_leading_imag_evl[i-1];
1801  }
1802 
1803  m_waveForceMag[0] = wavef_new;
1804 
1805  }
Nektar::VortexWaveInteraction::m_leading_imag_evl
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
Definition: VortexWaveInteraction.h:242
Vmath::Imin
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
Definition: Vmath.cpp:833
Nektar::VortexWaveInteraction::m_leading_real_evl
Array< OneD, NekDouble > m_leading_real_evl
Definition: VortexWaveInteraction.h:241
Nektar::VortexWaveInteraction::m_waveForceMag
Array< OneD, NekDouble > m_waveForceMag
Definition: VortexWaveInteraction.h:236
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047

Member Data Documentation

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_alpha
private

< Leading imaginary eigenvalue

Definition at line 244 of file VortexWaveInteraction.h.

Referenced by AppendEvlToFile(), CalcNonLinearWaveForce(), ExecuteLoop(), ExecuteWave(), UpdateAlpha(), UpdateDAlphaDWaveForceMag(), and VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_alphaStep
private

Definition at line 246 of file VortexWaveInteraction.h.

Referenced by UpdateAlpha(), and VortexWaveInteraction().

Array<OneD, Array<OneD, NekDouble > > Nektar::VortexWaveInteraction::m_bcsForcing
private

Definition at line 263 of file VortexWaveInteraction.h.

Referenced by FileRelaxation().

NekDouble Nektar::VortexWaveInteraction::m_dAlphaDWaveForceMag
private

Definition at line 250 of file VortexWaveInteraction.h.

Referenced by UpdateDAlphaDWaveForceMag(), and VortexWaveInteraction().

bool Nektar::VortexWaveInteraction::m_deltaFcnApprox
private

Definition at line 229 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_deltaFcnDecay
private

Definition at line 234 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_eigRelTol
private

Definition at line 248 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), and VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_iterEnd
private

Definition at line 222 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

bool Nektar::VortexWaveInteraction::m_iterinterface
private

Definition at line 253 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_iterStart
private

Definition at line 221 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_leading_imag_evl
private

< Leading real eigenvalue

Definition at line 242 of file VortexWaveInteraction.h.

Referenced by AppendEvlToFile(), CheckEigIsStationary(), CheckIfAtNeutralPoint(), ExecuteLoop(), ExecuteWave(), UpdateAlpha(), UpdateWaveForceMag(), and VortexWaveInteraction().

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_leading_real_evl
private

Definition at line 241 of file VortexWaveInteraction.h.

Referenced by AppendEvlToFile(), CheckEigIsStationary(), CheckIfAtNeutralPoint(), ExecuteWave(), UpdateAlpha(), UpdateWaveForceMag(), and VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_maxOuterIterations
private

Definition at line 225 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_maxWaveForceMagIter
private

Definition at line 227 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_minInnerIterations
private

Definition at line 226 of file VortexWaveInteraction.h.

Referenced by CheckEigIsStationary(), and VortexWaveInteraction().

bool Nektar::VortexWaveInteraction::m_moveMeshToCriticalLayer
private

Definition at line 232 of file VortexWaveInteraction.h.

Referenced by ExecuteLoop(), ExecuteRoll(), and VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_neutralPointTol
private

Definition at line 247 of file VortexWaveInteraction.h.

Referenced by CheckIfAtNeutralPoint(), and VortexWaveInteraction().

int Nektar::VortexWaveInteraction::m_nOuterIterations
private

Definition at line 224 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_prevAlpha
private

Definition at line 251 of file VortexWaveInteraction.h.

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_rollField
private

Definition at line 267 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and FileRelaxation().

NekDouble Nektar::VortexWaveInteraction::m_rollForceScale
private

Definition at line 239 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll(), and VortexWaveInteraction().

std::string Nektar::VortexWaveInteraction::m_sessionName
private

Definition at line 269 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), CopyFile(), ExecuteLoop(), ExecuteRoll(), FileRelaxation(), SaveLoopDetails(), and VortexWaveInteraction().

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionRoll
private

Definition at line 272 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), CheckIfAtNeutralPoint(), ExecuteLoop(), ExecuteRoll(), and VortexWaveInteraction().

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionStreak
private

Definition at line 275 of file VortexWaveInteraction.h.

Referenced by ExecuteStreak(), and VortexWaveInteraction().

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionVWI
private

Definition at line 270 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), CheckIfAtNeutralPoint(), ExecuteLoop(), ExecuteStreak(), ExecuteWave(), FileRelaxation(), GetReflectionIndex(), SaveLoopDetails(), VortexWaveInteraction(), and ~VortexWaveInteraction().

LibUtilities::SessionReaderSharedPtr Nektar::VortexWaveInteraction::m_sessionWave
private

Definition at line 276 of file VortexWaveInteraction.h.

Referenced by ExecuteWave(), and VortexWaveInteraction().

EquationSystemSharedPtr Nektar::VortexWaveInteraction::m_solverRoll
private

Definition at line 273 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), ExecuteRoll(), and VortexWaveInteraction().

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_streakField
private

Definition at line 265 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and ExecuteStreak().

bool Nektar::VortexWaveInteraction::m_useLinfPressureNorm
private

Definition at line 230 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), and VortexWaveInteraction().

Array<OneD, Array<OneD, NekDouble > > Nektar::VortexWaveInteraction::m_vwiForcing
private

Definition at line 260 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), ExecuteRoll(), and VortexWaveInteraction().

SolverUtils::ForcingProgrammaticSharedPtr Nektar::VortexWaveInteraction::m_vwiForcingObj
private

Definition at line 261 of file VortexWaveInteraction.h.

Referenced by ExecuteRoll().

VWIIterationType Nektar::VortexWaveInteraction::m_VWIIterationType
private

Definition at line 255 of file VortexWaveInteraction.h.

Referenced by VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_vwiRelaxation
private

Definition at line 249 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), FileRelaxation(), and VortexWaveInteraction().

Array<OneD, NekDouble> Nektar::VortexWaveInteraction::m_waveForceMag
private

Definition at line 236 of file VortexWaveInteraction.h.

Referenced by CalcNonLinearWaveForce(), UpdateWaveForceMag(), and VortexWaveInteraction().

NekDouble Nektar::VortexWaveInteraction::m_waveForceMagStep
private

Definition at line 237 of file VortexWaveInteraction.h.

Referenced by UpdateDAlphaDWaveForceMag(), UpdateWaveForceMag(), and VortexWaveInteraction().

MultiRegions::ExpListSharedPtr Nektar::VortexWaveInteraction::m_wavePressure
private

Definition at line 258 of file VortexWaveInteraction.h.

Referenced by CalcL2ToLinfPressure(), CalcNonLinearWaveForce(), and ExecuteWave().

Array<OneD, MultiRegions::ExpListSharedPtr> Nektar::VortexWaveInteraction::m_waveVelocities
private

Definition at line 257 of file VortexWaveInteraction.h.

Referenced by CalcL2ToLinfPressure(), CalcNonLinearWaveForce(), ExecuteWave(), and GetReflectionIndex().

Generated on Sat Jun 4 2016 00:02:35 for Nektar++ by   doxygen 1.8.8