47 VortexWaveInteraction::VortexWaveInteraction(
int argc,
char * argv[]):
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);
69 m_sessionVWI->LoadParameter(
"OuterIterationStoreSize",storesize,10);
115 if(
m_sessionVWI->DefinesSolverInfo(
"MoveMeshToCriticalLayer"))
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);
156 std::string vEquation =
m_sessionRoll->GetSolverInfo(
"SolverType");
166 int ncoeffs =
m_solverRoll->UpdateFields()[0]->GetNcoeffs();
169 for(
int i = 1; i < 4; ++i)
171 m_vwiForcing[i] = m_vwiForcing[i-1] + ncoeffs;
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);
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);
201 std::string LZstr(
"LZ");
203 cout <<
"Setting LZ in Linearised solver to " << LZ << endl;
209 std::string IterationTypeStr =
m_sessionVWI->GetSolverInfo(
"VWIIterationType");
228 m_sessionVWI->MatchSolverInfo(
"RestartIteration",
"True",restart,
false);
239 if((fp = fopen(
"OuterIter.his",
"r")))
242 std::vector<NekDouble> Alpha, Growth, Phase;
244 while(fgets(buf,BUFSIZ,fp))
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);
256 for(
int i = 0; i < nvals; ++i)
267 cout <<
" No File OuterIter.his to restart from" << endl;
273 string nstr = boost::lexical_cast<std::string>(
m_iterStart);
274 cout <<
"Restarting from iteration " <<
m_iterStart << endl;
276 cout <<
" " << rstfile << endl;
277 if(system(rstfile.c_str()))
281 std::string vwifile =
"cp -f Save/" + m_sessionName +
".vwi." + nstr +
" " + m_sessionName +
".vwi";
282 cout <<
" " << vwifile << endl;
283 if(system(vwifile.c_str()))
290 ASSERTL0(
false,
"Unknown VWIITerationType in restart");
297 if((fp = fopen(
"ConvergedSolns",
"r")))
300 std::vector<NekDouble> WaveForce, Alpha;
302 while(fgets(buf,BUFSIZ,fp))
304 sscanf(buf,
"%*d:%lf%lf",&waveforce,&alpha);
305 WaveForce.push_back(waveforce);
306 Alpha.push_back(alpha);
312 NekDouble min_alph = fabs(m_alpha[0]-Alpha[min_i]);
314 for(
int i = 1; i < Alpha.size(); ++i)
316 if(fabs(m_alpha[0]-Alpha[min_i]) < min_alph)
319 min_alph = fabs(m_alpha[0]-Alpha[min_i]);
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)
330 if(fabs(m_alpha[0]-Alpha[min_j]) < min_alph)
333 min_alph = fabs(m_alpha[0]-Alpha[min_j]);
338 if(fabs(Alpha[min_i] - Alpha[min_j]) > 1e-4)
358 string vEquation =
m_sessionRoll->GetSolverInfo(
"solvertype");
362 cout <<
"Executing Roll solver" << endl;
363 solverRoll->DoInitialise();
364 solverRoll->DoSolve();
365 solverRoll->Output();
367 for(
int g=0; g< solverRoll->GetNvariables(); ++g)
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;
379 string vEquation =
m_sessionRoll->GetSolverInfo(
"solvertype");
393 std::vector<std::string> vFieldNames =
m_sessionRoll->GetVariables();
394 vFieldNames.erase(vFieldNames.end()-1);
399 for(
int i = 0; i <
m_vwiForcingObj->UpdateForces().num_elements(); ++i)
413 for(
int i = 0; i <
m_vwiForcingObj->UpdateForces().num_elements(); ++i)
426 cout <<
"Executing Roll solver" << endl;
434 cout <<
"L 2 error (variable " <<
m_solverRoll->GetVariable(g) <<
") : " << vL2Error << endl;
435 cout <<
"L inf error (variable " <<
m_solverRoll->GetVariable(g) <<
") : " << vLinfError << endl;
442 cout <<
"Executing cp -f session.fld session.rst" << endl;
446 cout <<
"Writing data to session-Base.fld" << endl;
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();
455 outfield, variables);
463 std::string vDriverModule;
467 solverStreak->Execute();
475 cout <<
"Executing Streak Solver" << endl;
476 solverStreak->DoInitialise();
477 solverStreak->DoSolve();
478 solverStreak->Output();
484 for(
int g=0; g< solverStreak->GetNvariables(); ++g)
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;
494 cout <<
"Executing cp -f session.fld session_streak.fld" << endl;
502 std::string LZstr(
"LZ");
504 cout <<
"Setting LZ in Linearised solver to " << LZ << endl;
508 std::string vDriverModule;
509 m_sessionWave->LoadSolverInfo(
"Driver", vDriverModule,
"ModifiedArnoldi");
510 cout <<
"Setting up linearised NS sovler" << endl;
514 cout <<
"Executing wave solution " << endl;
515 solverWave->Execute();
518 cout <<
"Executing cp -f session_eig_0 session_eig_0.rst" << endl;
532 m_sessionWave->MatchSolverInfo(
"Driver",
"ModifiedArnoldi",defineshift,
true);
571 sprintf(c_alpha,
"%f",
m_alpha[0]);
576 syscall =
"../../utilities/PostProcessing/Extras/FldCalcBCs "
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()))
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()))
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()))
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()))
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()))
633 static int projectfield = -1;
644 if(projectfield == -1)
651 for(
int j = 0; j < BndConds.num_elements(); ++j)
659 if(projectfield != -1)
664 if(projectfield == -1)
666 cout <<
"using first field to project non-linear forcing which imposes a Dirichlet condition" << endl;
700 cout <<
"Area: " << area << endl;
701 invnorm = sqrt(area/invnorm);
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";
793 m_sessionVWI->MatchSolverInfo(
"Symmetrization",
"True",symm,
true);
816 for(i = 0; i <
npts; ++i)
818 coord[0] = coord_x[i];
819 coord[1] = coord_y[i];
827 for(i = 0; i <
npts; ++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);
845 for(i = 0; i <
npts; ++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);
859 m_waveVelocities[0]->GetPlane(0)->FwdTrans_BndConstrained(der2, m_vwiForceing[1]);
866 cout<<
"symmetrization is active"<<endl;
870 for(i = 0; i <
npts; ++i)
874 val[i] = 0.5*(der1[i] - der1[index[i]]);
884 for(i = 0; i <
npts; ++i)
888 val[i] = 0.5*(der2[i] - der2[index[i]]);
902 cout <<
"F_Linf: " <<
Vmath::Vmax(npts,val,1) << endl;
921 for(
int j = 0; j < 2; ++j)
925 for(
int i = 0; i <
npts; ++i)
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";
951 cout <<
"int P^2: " <<
m_wavePressure->GetPlane(0)->PhysIntegral(val) << endl;
953 cout <<
"PLinf: " <<
Vmath::Vmax(npts,val,1) << endl;
981 cout <<
"int P^2 " <<
m_wavePressure->GetPlane(0)->PhysIntegral(val) << endl;
986 cout <<
"Linf: " << Linf << endl;
998 l2 = sqrt(norm/area);
1000 cout <<
"L2: " << l2 << endl;
1002 cout <<
"Ratio Linf/L2: "<< Linf/l2 << endl;
1007 static map<string,int> opendir;
1009 if(opendir.find(dir) == opendir.end())
1012 string mkdir =
"mkdir " + dir;
1013 system(mkdir.c_str());
1018 string savefile = dir +
"/" + file +
"." + boost::lexical_cast<std::string>(n);
1019 string syscall =
"cp -f " + file +
" " + savefile;
1021 if(system(syscall.c_str()))
1031 static map<string,int> opendir;
1033 if(opendir.find(dir) == opendir.end())
1036 string mkdir =
"mkdir " + dir;
1037 system(mkdir.c_str());
1041 string savefile = dir +
"/" + file +
"." + boost::lexical_cast<std::string>(n);
1042 string syscall =
"mv -f " + file +
" " + savefile;
1044 if(system(syscall.c_str()))
1054 string syscall =
"cp -f " + cpfile1 +
" " + cpfile2;
1056 if(system(syscall.c_str()))
1065 fp = fopen(file.c_str(),
"a");
1073 fp = fopen(file.c_str(),
"a");
1107 bool skiprollstreak=
false;
1108 if(cnt==0 &&
m_sessionVWI->GetParameter(
"rollstreakfromit")==1)
1110 skiprollstreak =
true;
1111 cout<<
"skip roll-streak at the first iteration"<<endl;
1115 if(skiprollstreak !=
true)
1133 string movedinterpmesh =
m_sessionName +
"_interp_moved.xml";
1143 sprintf(c,
"%d",cnt);
1146 syscall =
"cp -f " +
m_sessionName+
"-Base.fld" +
" " + oldroll;
1147 cout<<syscall.c_str()<<endl;
1148 if(system(syscall.c_str()))
1156 string filewavepressure =
m_sessionName +
"_wave_p_split.fld";
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]);
1165 if(
m_sessionVWI->GetSolverInfo(
"INTERFACE")!=
"phase" )
1167 cout<<
"zerophase"<<endl;
1169 syscall =
"../../utilities/PostProcessing/Extras/MoveMesh "
1170 + filePost +
" "+ filestreak +
" "+ fileinterp +
" "+ alpchar;
1172 cout<<syscall.c_str()<<endl;
1173 if(system(syscall.c_str()))
1179 syscall =
"../../utilities/PostProcessing/Extras/MoveMesh "
1180 + filePost +
" " + filestreak +
" " + filePost +
" "+ alpchar;
1181 cout<<syscall.c_str()<<endl;
1182 if(system(syscall.c_str()))
1191 syscall =
"cp -f " + filestreak +
" " + oldstreak;
1192 cout<<syscall.c_str()<<endl;
1193 if(system(syscall.c_str()))
1200 string movedinterpmesh =
m_sessionName +
"_interp_moved.xml";
1204 syscall =
"../../utilities/PostProcessing/Extras/FieldToField "
1205 + fileinterp +
" " + filestreak +
" " + movedinterpmesh +
" "
1208 cout<<syscall.c_str()<<endl;
1209 if(system(syscall.c_str()))
1218 syscall =
"cp -f " + meshfile +
" " + meshold;
1219 cout<<syscall.c_str()<<endl;
1220 if(system(syscall.c_str()))
1226 syscall =
"cp -f " + movedmesh +
" " + meshfile;
1227 cout<<syscall.c_str()<<endl;
1228 if(system(syscall.c_str()))
1234 syscall =
"cp -f " + interpstreak +
" " + filestreak;
1235 cout<<syscall.c_str()<<endl;
1236 if(system(syscall.c_str()))
1247 cout<<syscall.c_str()<<endl;
1248 if(system(syscall.c_str()))
1255 syscall =
"cp -f " + ujump +
" " +
m_sessionName+
"_u_5.bc_"+c;
1256 cout<<syscall.c_str()<<endl;
1257 if(system(syscall.c_str()))
1263 syscall =
"cp -f " + vjump +
" " +
m_sessionName+
"_v_5.bc_"+c;
1264 cout<<syscall.c_str()<<endl;
1265 if(system(syscall.c_str()))
1283 sprintf(c1,
"%d",cnt);
1286 syscall =
"../../utilities/PostProcessing/Extras/FldCalcBCs "
1287 + movedmesh +
" " + wavefile +
" " + interpstreak +
"> data"+c1;
1288 cout<<syscall.c_str()<<endl;
1289 if(system(syscall.c_str()))
1295 syscall =
"cp -f " + movedinterpmesh +
" " + fileinterp;
1296 cout<<syscall.c_str()<<endl;
1297 if(system(syscall.c_str()))
1302 syscall =
"cp -f " + movedmesh +
" " + filePost;
1303 cout<<syscall.c_str()<<endl;
1304 if(system(syscall.c_str()))
1311 else if(
m_sessionVWI->GetSolverInfo(
"INTERFACE")==
"phase" )
1313 cout<<
"phase"<<endl;
1324 syscall =
"cp -f " + filestreak +
" " + oldstreak;
1325 cout<<syscall.c_str()<<endl;
1326 if(system(syscall.c_str()))
1333 cout<<syscall.c_str()<<endl;
1334 if(system(syscall.c_str()))
1341 syscall =
"cp -f " + meshfile +
" " + meshold;
1342 cout<<syscall.c_str()<<endl;
1343 if(system(syscall.c_str()))
1351 cout<<syscall.c_str()<<endl;
1352 if(system(syscall.c_str()))
1359 syscall =
"cp -f " + ujump +
" " +
m_sessionName+
"_u_5.bc_"+c;
1360 cout<<syscall.c_str()<<endl;
1361 if(system(syscall.c_str()))
1367 syscall =
"cp -f " + vjump +
" " +
m_sessionName+
"_v_5.bc_"+c;
1368 cout<<syscall.c_str()<<endl;
1369 if(system(syscall.c_str()))
1379 cout<<
"cr="<<cr_str<<endl;
1382 syscall =
"../../utilities/PostProcessing/Extras/MoveMesh "
1383 + filePost +
" "+ filestreak +
" "+ fileinterp +
" "+ alpchar
1386 cout<<syscall.c_str()<<endl;
1387 if(system(syscall.c_str()))
1393 syscall =
"../../utilities/PostProcessing/Extras/MoveMesh "
1394 + filePost +
" " + filestreak +
" " + filePost +
" "+ alpchar
1396 cout<<syscall.c_str()<<endl;
1397 if(system(syscall.c_str()))
1403 syscall =
"../../utilities/PostProcessing/Extras/FieldToField "
1404 + fileinterp +
" " + filestreak +
" " + movedinterpmesh +
" "
1407 cout<<syscall.c_str()<<endl;
1408 if(system(syscall.c_str()))
1414 syscall =
"../../utilities/PostProcessing/Extras/SplitFld "
1415 + filePost +
" " + filewave;
1417 cout<<syscall.c_str()<<endl;
1418 if(system(syscall.c_str()))
1423 syscall =
"../../utilities/PostProcessing/Extras/FieldToField "
1424 + filePost +
" " + filewavepressure +
" " + movedmesh +
" "
1425 + interwavepressure;
1427 cout<<syscall.c_str()<<endl;
1428 if(system(syscall.c_str()))
1445 sprintf(c1,
"%d",cnt);
1449 syscall =
"cp -f "+ interwavepressure +
" "+
m_sessionName+
".fld";
1450 cout<<syscall.c_str()<<endl;
1451 if(system(syscall.c_str()))
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()))
1469 syscall =
"cp -f " + movedmesh +
" " + meshfile;
1470 cout<<syscall.c_str()<<endl;
1471 if(system(syscall.c_str()))
1477 syscall =
"cp -f " + interpstreak +
" " + filestreak;
1478 cout<<syscall.c_str()<<endl;
1479 if(system(syscall.c_str()))
1484 syscall =
"cp -f " + movedinterpmesh +
" " + fileinterp;
1485 cout<<syscall.c_str()<<endl;
1486 if(system(syscall.c_str()))
1491 syscall =
"cp -f " + movedmesh +
" " + filePost;
1492 cout<<syscall.c_str()<<endl;
1493 if(system(syscall.c_str()))
1516 char alpchar[16]=
"";
1517 sprintf(alpchar,
"%f",
m_alpha[0]);
1522 string filewavepressure =
m_sessionName +
"_wave_p_split.fld";
1525 string interwavepressure =
m_sessionName +
"_wave_p_split_interp.fld";
1526 syscall =
"../../utilities/PostProcessing/Extras/MoveMesh "
1527 + filePost +
" "+ filestreak +
" "+ fileinterp +
" "+ alpchar;
1529 cout<<syscall.c_str()<<endl;
1530 if(system(syscall.c_str()))
1536 syscall =
"../../utilities/PostProcessing/Extras/MoveMesh "
1537 + filePost +
" " + filestreak +
" " + filePost +
" "+ alpchar;
1538 cout<<syscall.c_str()<<endl;
1539 if(system(syscall.c_str()))
1546 syscall =
"cp -f " + filestreak +
" " + oldstreak;
1547 cout<<syscall.c_str()<<endl;
1548 if(system(syscall.c_str()))
1555 string movedinterpmesh =
m_sessionName +
"_interp_moved.xml";
1559 syscall =
"../../utilities/PostProcessing/Extras/FieldToField "
1560 + fileinterp +
" " + filestreak +
" " + movedinterpmesh
1561 +
" " + interpstreak;
1563 cout<<syscall.c_str()<<endl;
1564 if(system(syscall.c_str()))
1572 syscall =
"cp -f " + meshfile +
" " + meshold;
1573 cout<<syscall.c_str()<<endl;
1574 if(system(syscall.c_str()))
1580 syscall =
"cp -f " + movedmesh +
" " + meshfile;
1581 cout<<syscall.c_str()<<endl;
1582 if(system(syscall.c_str()))
1588 syscall =
"cp -f " + interpstreak +
" " + filestreak;
1589 cout<<syscall.c_str()<<endl;
1590 if(system(syscall.c_str()))
1607 static NekDouble previous_real_evl = -1.0;
1608 static NekDouble previous_imag_evl = -1.0;
1609 static int min_iter = 0;
1613 previous_real_evl = -1.0;
1642 cout <<
"Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1655 bool returnval =
false;
1685 cout <<
"Warning: imaginary eigenvalue is greater than 1e-2" << endl;
1723 for(i = 0; i < nstore; ++i)
1727 store = WaveForce[i];
1728 WaveForce[i] = WaveForce[i+k];
1729 WaveForce[i+k] = store;
1732 Growth[i] = Growth[i+k];
1733 Growth[i+k] = store;
1737 for(i = 0; i < nstore-1; ++i)
1739 if(Growth[i]*Growth[i+1] < 0.0)
1749 wavef_new = (WaveForce[0]*Growth[1] - WaveForce[1]*Growth[0])/(Growth[1]-Growth[0]);
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]);
1767 for(j = 1; j < nsteps+1; ++j)
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]);
1774 if(gval*gval_m1 < 0.0)
1776 wavef_new = ((a+da/(
NekDouble)nsteps)*gval - a*gval_m1)/
1832 int nstore = (
m_alpha.num_elements() < outeriter)?
m_alpha.num_elements(): outeriter;
1842 for(i = 0; i < nstore; ++i)
1847 Alpha[i] = Alpha[i+k];
1851 Growth[i] = Growth[i+k];
1852 Growth[i+k] = store;
1856 for(i = 0; i < nstore-1; ++i)
1858 if(Growth[i]*Growth[i+1] < 0.0)
1868 alp_new = (Alpha[0]*Growth[1] - Alpha[1]*Growth[0])/(Growth[1]-Growth[0]);
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]);
1886 for(j = 1; j < nsteps+1; ++j)
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]);
1893 if(gval*gval_m1 < 0.0)
1895 alp_new = ((a+da/(
NekDouble)nsteps)*gval - a*gval_m1)/
1915 for(
int i =
m_alpha.num_elements()-1; i > 0; --i)
1946 bool useOnlyQuads =
false;
1947 if(
m_sessionVWI->DefinesSolverInfo(
"SymmetriseOnlyQuads"))
1949 useOnlyQuads =
true;
1953 for(
int e = 0; e < nel; ++e)
1962 for(i = 0; i < e_npts; ++i)
1970 for(i = cnt; i < cnt+e_npts; ++i)
1972 xnew = - coord_x[i] + xmax;
1973 ynew = - coord_y[i];
1975 for(j = start; j >=0 ; --j)
1977 if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1988 for(j = npts-1; j > start; --j)
1991 if((coord_x[j]-xnew)*(coord_x[j]-xnew) + (coord_y[j]-ynew)*(coord_y[j]-ynew) < tol)
1997 ASSERTL0(j != start,
"Failed to find matching point");
2007 cout<<
"relaxation..."<<endl;
2015 *boost::static_pointer_cast<MultiRegions::ExpList1D>(Iexp[reg]));
2016 int nq = Ilayer->GetTotPoints();
2021 for(
int i = 1; i < 4; ++i)
2023 m_bcsForcing[i] = m_bcsForcing[i-1] + nq;
2036 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef_u;
2037 std::vector<std::vector<NekDouble> > FieldData_u;
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());
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);
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());
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);
static void DisableManagement(std::string whichPool="")
#define ASSERTL0(condition, msg)
LibUtilities::SessionReaderSharedPtr m_sessionRoll
MultiRegions::ExpListSharedPtr m_wavePressure
void CalcL2ToLinfPressure(void)
boost::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
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.
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
void ExecuteLoop(bool CalcWaveForce=true)
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
std::string m_sessionName
Array< OneD, NekDouble > m_leading_imag_evl
< Leading real eigenvalue
int Imin(int n, const T *x, const int incx)
Return the index of the minimum element in x.
boost::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
Array< OneD, MultiRegions::ExpListSharedPtr > m_rollField
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
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
ForcingFactory & GetForcingFactory()
Declaration of the forcing factory singleton.
bool CheckEigIsStationary(bool reset=false)
VWIIterationType GetVWIIterationType(void)
bool m_moveMeshToCriticalLayer
void SaveFile(std::string fileend, std::string dir, int n)
Array< OneD, MultiRegions::ExpListSharedPtr > m_waveVelocities
void UpdateDAlphaDWaveForceMag(NekDouble alphainit)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
void MoveFile(std::string fileend, std::string dir, int n)
EquationSystemSharedPtr m_solverRoll
SolverUtils::ForcingProgrammaticSharedPtr m_vwiForcingObj
void UpdateWaveForceMag(int n)
static void EnableManagement(std::string whichPool="")
boost::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
Array< OneD, NekDouble > m_leading_real_evl
NekDouble m_vwiRelaxation
NekDouble m_rollForceScale
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
Array< OneD, Array< OneD, NekDouble > > m_bcsForcing
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
NekDouble m_neutralPointTol
LibUtilities::SessionReaderSharedPtr m_sessionWave
~VortexWaveInteraction(void)
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
const std::string VWIIterationTypeMap[]
Array< OneD, NekDouble > m_waveForceMag
VWIIterationType m_VWIIterationType
void SaveLoopDetails(std::string dir, int i)
void Neg(int n, T *x, const int incx)
Negate x = -x.
boost::shared_ptr< FieldIO > FieldIOSharedPtr
bool m_useLinfPressureNorm
void FileRelaxation(int reg)
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
NekDouble m_waveForceMagStep
EquationSystemFactory & GetEquationSystemFactory()
Array< OneD, Array< OneD, NekDouble > > m_vwiForcing
LibUtilities::SessionReaderSharedPtr m_sessionVWI
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.
void CalcNonLinearWaveForce(void)
NekDouble m_dAlphaDWaveForceMag
void AddForcing(const SolverUtils::ForcingSharedPtr &pForce)
int m_maxWaveForceMagIter
void CopyFile(std::string file1end, std::string file2end)
This class is the base class for Navier Stokes problems.
DriverFactory & GetDriverFactory()
NekDouble m_deltaFcnDecay
Array< OneD, MultiRegions::ExpListSharedPtr > m_streakField
Array< OneD, NekDouble > m_alpha
< Leading imaginary eigenvalue
Array< OneD, int > GetReflectionIndex(void)
void AppendEvlToFile(std::string file, int n)
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
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.
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.
bool CheckIfAtNeutralPoint(void)
LibUtilities::SessionReaderSharedPtr m_sessionStreak