Commit 030d2299 authored by Andrea Isoni's avatar Andrea Isoni

update utilities and VWI


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@3668 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent ef44d7f3
SET(VortexWaveInteractionSolverSource
# ../Auxiliary/Driver.cpp
# ../Auxiliary/DriverStandard.cpp
# ../Auxiliary/DriverArnoldi.cpp
# ../Auxiliary/DriverModifiedArnoldi.cpp
# ../Auxiliary/EquationSystem.cpp
# ../ADRSolver/EquationSystems/SteadyAdvectionDiffusion.cpp
# ../ADRSolver/EquationSystems/SteadyAdvectionDiffusionReaction.cpp
# ../IncNavierStokesSolver/EquationSystems/CoupledLinearNS.cpp
# ../IncNavierStokesSolver/EquationSystems/CoupledLocalToGlobalC0ContMap.cpp
# ../IncNavierStokesSolver/EquationSystems/IncNavierStokes.cpp
# ../IncNavierStokesSolver/EquationSystems/VelocityCorrectionScheme.cpp
# ../IncNavierStokesSolver/AdvectionTerms/AdvectionTerm.cpp
# ../IncNavierStokesSolver/AdvectionTerms/NavierStokesAdvection.cpp
# ../IncNavierStokesSolver/AdvectionTerms/LinearisedAdvection.cpp
# ../IncNavierStokesSolver/AdvectionTerms/AdjointAdvection.cpp
# ./VortexWaveInteraction.cpp
# ./VortexWaveInteractionSolver.cpp
../Auxiliary/Driver.cpp
../Auxiliary/DriverStandard.cpp
../Auxiliary/DriverArnoldi.cpp
../Auxiliary/DriverModifiedArnoldi.cpp
../Auxiliary/EquationSystem.cpp
../ADRSolver/EquationSystems/SteadyAdvectionDiffusion.cpp
../ADRSolver/EquationSystems/SteadyAdvectionDiffusionReaction.cpp
../IncNavierStokesSolver/EquationSystems/CoupledLinearNS.cpp
../IncNavierStokesSolver/EquationSystems/CoupledLocalToGlobalC0ContMap.cpp
../IncNavierStokesSolver/EquationSystems/IncNavierStokes.cpp
../IncNavierStokesSolver/EquationSystems/VelocityCorrectionScheme.cpp
../IncNavierStokesSolver/AdvectionTerms/AdvectionTerm.cpp
../IncNavierStokesSolver/AdvectionTerms/NavierStokesAdvection.cpp
../IncNavierStokesSolver/AdvectionTerms/LinearisedAdvection.cpp
../IncNavierStokesSolver/AdvectionTerms/AdjointAdvection.cpp
./VortexWaveInteraction.cpp
./VortexWaveInteractionSolver.cpp
)
IF (NEKTAR_USE_ARPACK)
......@@ -23,6 +23,6 @@ IF (NEKTAR_USE_ARPACK)
../Auxiliary/DriverArpack.cpp)
ENDIF (NEKTAR_USE_ARPACK)
#ADD_SOLVER_EXECUTABLE(VortexWaveInteractionSolver solvers
# ${VortexWaveInteractionSolverSource})
ADD_SOLVER_EXECUTABLE(VortexWaveInteractionSolver solvers
${VortexWaveInteractionSolverSource})
......@@ -39,7 +39,7 @@ int main(int argc, char *argv[])
MultiRegions::ContField1DSharedPtr &outfieldx,
MultiRegions::ContField1DSharedPtr &outfieldy,
MultiRegions::ExpListSharedPtr &streak,
Array<OneD, int> Refindices);
Array<OneD, int> Refindices, NekDouble alpha);
void Manipulate(Array<OneD, int> Iregions, int coordim, SpatialDomains::MeshGraphSharedPtr &mesh,
SpatialDomains::BoundaryConditions &bcs,
Array<OneD,MultiRegions::ExpList1DSharedPtr> &infields,
......@@ -61,17 +61,20 @@ decomment the lines which follow '//decomment'
int i,j;
if(argc != 4)
if(argc >= 6 || argc <4)
{
fprintf(stderr,"Usage: ./FldCalcBCs meshfile fieldfile streakfile\n");
fprintf(stderr,"Usage: ./FldCalcBCs meshfile fieldfile streakfile (alpha)\n");
exit(1);
}
cout<<"argc="<<argc<<endl;
//we consider by default the case when alpha is provided: argc=5
if(argc==4){ argc=5;}
//----------------------------------------------
string meshfile(argv[argc-3]);
string meshfile(argv[argc-4]);
//setting argc=2 will take consider only the first two words argv[0],argv[1]
LibUtilities::SessionReaderSharedPtr vSession
= LibUtilities::SessionReader::CreateInstance(argc, argv);
= LibUtilities::SessionReader::CreateInstance(2, argv);
// Read in mesh from input file
SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(meshfile);
......@@ -84,9 +87,20 @@ decomment the lines which follow '//decomment'
::AllocateSharedPtr(vSession, graphShPt);
SpatialDomains::BoundaryConditions bcs(vSession, graphShPt);
//----------------------------------------------
//load alpha is provided
NekDouble alpha=0;
string alp_s;
if(argv[argc-1])
{
alp_s = argv[argc-1];
alpha = boost::lexical_cast<double>(alp_s);
}
//----------------------------------------------
// Import field file.
string fieldfile(argv[argc-2]);
string fieldfile(argv[argc-3]);
vector<SpatialDomains::FieldDefinitionsSharedPtr> fielddef;
vector<vector<NekDouble> > fielddata;
graphShPt->Import(fieldfile,fielddef,fielddata);
......@@ -96,7 +110,7 @@ decomment the lines which follow '//decomment'
int nfields;
nfields = fielddef[0]->m_fields.size();
Array<OneD, MultiRegions::ExpListSharedPtr> fields;
fields= Array<OneD, MultiRegions::ExpListSharedPtr>(nfields);
fields= Array<OneD, MultiRegions::ExpListSharedPtr>(nfields);
std::string solvtype = vSession->GetSolverInfo("SOLVERTYPE");
......@@ -175,7 +189,7 @@ decomment the lines which follow '//decomment'
//streak = MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr
// (vSession, graphShPt, true, "w");
string streakfile(argv[argc-1]);
string streakfile(argv[argc-2]);
vector<SpatialDomains::FieldDefinitionsSharedPtr> streakdef;
vector<vector<NekDouble> > streakdata;
graphShPt->Import(streakfile,streakdef,streakdata);
......@@ -246,7 +260,7 @@ cout<<"OOOK"<<endl;
int coordim = graphShPt->GetMeshDimension();
//remark Ilayers[2] is the critical layer
static Array<OneD, int> Refindices = GetReflectionIndex(fields, Ilayers[2]);
Extractlayerdata(Ilayers,coordim, graphShPt,vSession, bcs, fields, outfieldx,outfieldy,streak,Refindices);
Extractlayerdata(Ilayers,coordim, graphShPt,vSession, bcs, fields, outfieldx,outfieldy,streak,Refindices, alpha);
//--------------------------------------------------------------------------------------
......@@ -447,10 +461,10 @@ cout<<"OOOK"<<endl;
NekDouble xnew,ynew;
int start = npts-1;
cout<<"xmax="<<xmax<<endl;
//cout<<"xmax="<<xmax<<endl;
for(i = 0; i < npts; ++i)
{
cout<<"Ref index for x="<<coord_x[i]<<endl;
//cout<<"Ref index for x="<<coord_x[i]<<endl;
xnew = - coord_x[i] + xmax;
ynew = - coord_y[i];
......@@ -490,7 +504,7 @@ cout<<"Ref index for x="<<coord_x[i]<<endl;
MultiRegions::ContField1DSharedPtr &outfieldx,
MultiRegions::ContField1DSharedPtr &outfieldy,
MultiRegions::ExpListSharedPtr &streak,
Array<OneD, int> Refindices)
Array<OneD, int> Refindices, NekDouble alpha)
{
//1 I regions is expected: (the layer is the last region)
ASSERTL0(Iregions.num_elements()==3, "something wrong with the number of I layers");
......@@ -1161,8 +1175,11 @@ cout<<"x"<<" P_re"<<" dP_re"<<" streak"<<" dstreak"<<" pjump"<<endl;
NekDouble rho ;
//rho = 0.08;
rho = session->GetParameter("RHO");
NekDouble alpha;
alpha = session->GetParameter("ALPHA");
//load alpha only if not manually specified
if(alpha==0)
{
alpha = session->GetParameter("ALPHA");
}
cout<<"alpha="<<alpha<<endl;
NekDouble alpha53;
......@@ -1316,7 +1333,7 @@ cout<<"symmetrise the jump conditions"<<endl;
NekDouble jactest;
NekDouble laytest=0;
for(int g=0; g<nq1D; g++)
{
//pjump[g] = n0*curv[g]*mu53[g]*dP_square[g] ;
......@@ -1346,7 +1363,7 @@ outfieldx->GetPhys()[g]<<" "
<<outfieldy->GetPhys()[g]
<<endl;
*/
laytest += stphysreg[g];
cout<<setw(14)<<x0d[g]<<" "<<setw(13)<<x1d[g]<<
" "<<Rephysreg[g]<<" "<<dP_re[g]<<" "
......@@ -1356,7 +1373,7 @@ cout<<setw(14)<<x0d[g]<<" "<<setw(13)<<x1d[g]<<
//(gmat0[g]*tx[g] +gmat2[g]*ty[g])<<" "<<1/Jac[g]
//<<" "<<dermu53[g]
//<<" "<<derfun1D[g]
tx[g]<<" "<<ty[g]<<
tx[g]<<" "<<ty[g]<<" "<<
//f1[g]<<" "<<f2[g]
//prod_b[g]<<" "<<dersprod_b[g]
//fun1D[g]<<" "<<derfun1D[g]<<" "<<
......@@ -1366,10 +1383,11 @@ pjump[g]<<setw(13)<<" "<<d2v[g]<<" "
<<outfieldx->GetPhys()[g]<<" "
<<outfieldy->GetPhys()[g]<<endl;
//cout<<"jactest="<<jactest<<endl;
ASSERTL0(abs(jactest)<0.5, "jac 1D problem..");
ASSERTL0(abs(jactest)<0.5, "jac 1D problem..");
}
ASSERTL0(abs(laytest)<0.001, "critical layer wrong");
// gamma(1/3)= 2.6789385347077476337
// need the tangents related to the expList1D outfieldx[region]
......
......@@ -58,14 +58,14 @@ int main(int argc, char *argv[])
Array<OneD, NekDouble> x_lay, Array<OneD, NekDouble> y_lay,
Array<OneD, NekDouble> & Pcurvx,
Array<OneD, NekDouble> & Pcurvy,
Array<OneD, int>&Eids, int Npoints);
Array<OneD, int>&Eids, int Npoints, string s_alp);
int i,j;
if(argc != 4)
if(argc != 5)
{
fprintf(stderr,"Usage: ./MoveMesh meshfile fieldfile changefile\n");
fprintf(stderr,"Usage: ./MoveMesh meshfile fieldfile changefile alpha\n");
exit(1);
}
//ATTEnTION !!! with argc=2 you impose that vSession refers to is argv[1]=meshfile!!!!!
......@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
//----------------------------------------------
// Read in mesh from input file
string meshfile(argv[argc-3]);
string meshfile(argv[argc-4]);
SpatialDomains::MeshGraphSharedPtr graphShPt = SpatialDomains::MeshGraph::Read(meshfile);
//----------------------------------------------
......@@ -98,15 +98,23 @@ int main(int argc, char *argv[])
// store name of the file to change
string changefile(argv[argc-1]);
string changefile(argv[argc-2]);
//----------------------------------------------
//store the value of alpha
string charalp (argv[argc-1]);
//NekDouble alpha = boost::lexical_cast<double>(charalp);
cout<<"read alpha="<<charalp<<endl;
//----------------------------------------
// Import field file.
string fieldfile(argv[argc-2]);
string fieldfile(argv[argc-3]);
vector<SpatialDomains::FieldDefinitionsSharedPtr> fielddef;
vector<vector<NekDouble> > fielddata;
graphShPt->Import(fieldfile,fielddef,fielddata);
//----------------------------------------------
//cout<<"dim st="<<fielddef[0]->m_fields.size()<<endl;
//fill a vector with the streak sol
......@@ -850,7 +858,7 @@ cout<<i<<" "<<xnew[i]<<" "<<ynew[i]<<endl;
//replace the vertices with the new ones
//Replacevertices(changefile, xnew , ynew, x_c, y_c, Cpointsx, Cpointsy, Eids, npedge);
Replacevertices(changefile, xnew , ynew, x_c, y_c, Addpointsx, Addpointsy, Eids, npedge);
Replacevertices(changefile, xnew , ynew, x_c, y_c, Addpointsx, Addpointsy, Eids, npedge,charalp);
}
......@@ -1568,7 +1576,7 @@ cout<<"gen newton xi="<<xi<<" yi="<<yi<<" elmtid="<<elmtid<<endl;
Array<OneD, NekDouble> x_lay, Array<OneD, NekDouble> y_lay,
Array<OneD, NekDouble> & Pcurvx,
Array<OneD, NekDouble> & Pcurvy,
Array<OneD, int>&Eids, int Npoints)
Array<OneD, int>&Eids, int Npoints, string s_alp)
{
//load existing file
string newfile;
......@@ -1605,12 +1613,62 @@ cout<<"gen newton xi="<<xi<<" yi="<<yi<<" elmtid="<<elmtid<<endl;
TiXmlHandle docHandlenew(&docnew);
TiXmlNode* nodenew = NULL;
TiXmlElement* meshnew = NULL;
TiXmlElement* masternew = NULL; // Master tag within which all data is contained.
TiXmlElement* masternew = NULL;
TiXmlElement* condnew = NULL;
TiXmlElement* Parsnew = NULL;
TiXmlElement* parnew = NULL;
// Master tag within which all data is contained.
masternew = docnew.FirstChildElement("NEKTAR");
ASSERTL0(masternew, "Unable to find NEKTAR tag in file.");
//set the alpha value
string alphastring;
condnew = masternew->FirstChildElement("CONDITIONS");
Parsnew = condnew->FirstChildElement("PARAMETERS");
cout<<"alpha="<<s_alp<<endl;
parnew = Parsnew->FirstChildElement("P");
while(parnew)
{
TiXmlNode *node = parnew->FirstChild();
if (node)
{
// Format is "paramName = value"
std::string line = node->ToText()->Value();
std::string lhs;
std::string rhs;
/// Pull out lhs and rhs and eliminate any spaces.
int beg = line.find_first_not_of(" ");
int end = line.find_first_of("=");
// Check for no parameter name
if (beg == end) throw 1;
// Check for no parameter value
if (end != line.find_last_of("=")) throw 1;
// Check for no equals sign
if (end == std::string::npos) throw 1;
lhs = line.substr(line.find_first_not_of(" "), end-beg);
lhs = lhs.substr(0, lhs.find_last_not_of(" ")+1);
//rhs = line.substr(line.find_last_of("=")+1);
//rhs = rhs.substr(rhs.find_first_not_of(" "));
//rhs = rhs.substr(0, rhs.find_last_not_of(" ")+1);
boost::to_upper(lhs);
if(lhs == "ALPHA")
{
alphastring = "Alpha = "+s_alp;
parnew->RemoveChild(node);
parnew->LinkEndChild(new TiXmlText(alphastring));
}
}
parnew = parnew->NextSiblingElement("P");
}
// Find the Mesh tag and same the dim and space attributes
meshnew = masternew->FirstChildElement("GEOMETRY");
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment