Commit 56d14b50 authored by Sergey Yakovlev's avatar Sergey Yakovlev

Changed TimingCGHelmSolve2D

Also added SetSolveInfo to SessionReader
Need to change HDG timing in 2D and adjust scripts accordingly
parent c00066a1
......@@ -57,9 +57,11 @@ int main(int argc, char *argv[])
if (vComm->GetRank() == 0)
{
cout << "Solving 2D Helmholtz:" << endl;
cout << " Lambda : " << factors[StdRegions::eFactorLambda] << endl;
cout << " No. modes : " << bkey0.GetNumModes() << endl;
cout << "Solving 2D Helmholtz: " << endl;
cout << " Communication: " << vSession->GetComm()->GetType() << endl;
cout << " Solver type : " << vSession->GetSolverInfo("GlobalSysSoln") << endl;
cout << " Lambda : " << factors[StdRegions::eFactorLambda] << endl;
cout << " No. modes : " << bkey0.GetNumModes() << endl;
cout << endl;
}
......
......@@ -57,10 +57,18 @@ int main(int argc, char *argv[])
if (vComm->GetRank() == 0)
{
cout << "Solving 3D Helmholtz:" << endl;
cout << " Lambda : " << factors[StdRegions::eFactorLambda] << endl;
cout << " No. modes : " << bkey0.GetNumModes() << endl;
cout << endl;
cout << "Solving 3D Helmholtz:" << endl;
cout << " - Communication: "
<< vSession->GetComm()->GetType() << " ("
<< vSession->GetComm()->GetSize()
<< " processes)" << endl;
cout << " - Solver type : "
<< vSession->GetSolverInfo("GlobalSysSoln") << endl;
cout << " - Lambda : "
<< factors[StdRegions::eFactorLambda] << endl;
cout << " - No. modes : "
<< bkey0.GetNumModes() << endl;
cout << endl;
}
//----------------------------------------------
......
......@@ -669,6 +669,20 @@ namespace Nektar
return iter->second;
}
/**
*
*/
void SessionReader::SetSolverInfo(
const std::string &pProperty, const std::string &pValue)
{
std::string vProperty = boost::to_upper_copy(pProperty);
SolverInfoMap::iterator iter = m_solverInfo.find(vProperty);
ASSERTL1(iter != m_solverInfo.end(),
"Unable to find requested property: " + pProperty);
iter->second = pValue;
}
/**
*
......
......@@ -227,6 +227,9 @@ namespace Nektar
/// Returns the value of the specified solver info property.
LIB_UTILITIES_EXPORT const std::string& GetSolverInfo(
const std::string &pProperty) const;
/// Sets the value of the specified solver info property.
LIB_UTILITIES_EXPORT void SetSolverInfo(
const std::string &pProperty, const std::string &pValue);
/// Returns the value of the specified solver info property as enum
template<typename T>
inline const T GetSolverInfoAsEnum(const std::string &pName) const;
......
......@@ -86,7 +86,7 @@ namespace Nektar
const std::string &variable,
const bool SetUpJustDG,
const bool DeclareCoeffPhysArrays)
: ExpList2D(pSession, graph2D, DeclareCoeffPhysArrays),
: ExpList2D(pSession,graph2D,DeclareCoeffPhysArrays,variable),
m_bndCondExpansions(),
m_bndConditions(),
m_trace(NullExpListSharedPtr),
......
......@@ -3,6 +3,7 @@
<CONDITIONS>
<PARAMETERS>
<P> w = 3 </P>
<P> Lambda = 1 </P>
</PARAMETERS>
......@@ -12,15 +13,22 @@
<!--These composites must be defined in the geometry file.-->
<BOUNDARYREGIONS>
<B ID="0"> C[1-4] </B>
<!--
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
<B ID="2"> C[3] </B>
<B ID="3"> C[4] </B>
<B ID="3"> C[4] </B>
-->
</BOUNDARYREGIONS>
<!--The region numbers below correspond to the regions specified in the
BoundaryRegion definition above.-->
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" VALUE="sin(w*PI*x)*sin(w*PI*y)" />
</REGION>
<!--
<REGION REF="0">
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
</REGION>
......@@ -33,14 +41,15 @@
<REGION REF="3">
<D VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
</REGION>
-->
</BOUNDARYCONDITIONS>
<FUNCTION NAME="Forcing">
<E VAR="u" VALUE="-(Lambda + 2*PI*PI)*sin(PI*x)*sin(PI*y)" />
<E VAR="u" VALUE="-(Lambda + 2*w*w*PI*PI)*sin(PI*x)*sin(PI*y)" />
</FUNCTION>
<FUNCTION NAME="ExactSolution">
<E VAR="u" VALUE="sin(PI*x)*sin(PI*y)" />
<E VAR="u" VALUE="sin(w*PI*x)*sin(w*PI*y)" />
</FUNCTION>
......
......@@ -18,389 +18,392 @@ std::string PortablePath(const boost::filesystem::path& path);
int main(int argc, char *argv[])
{
MultiRegions::ContField2DSharedPtr Exp,Fce,Sol;
int i, nq, coordim;
Array<OneD,NekDouble> fce,sol;
Array<OneD,NekDouble> xc0,xc1,xc2;
NekDouble lambda;
vector<string> vFilenames;
if(argc != 5)
{
fprintf(stderr,"Usage: TimingCGHelmSolve2D Type MeshSize NumModes OptimisationLevel\n");
fprintf(stderr," where: - Type is one of the following:\n");
fprintf(stderr," 1: Regular Quadrilaterals \n");
fprintf(stderr," 2: Deformed Quadrilaterals (may not be supported) \n");
fprintf(stderr," 3: Regular Triangles \n");
fprintf(stderr," where: - MeshSize is 1/h \n");
fprintf(stderr," where: - NumModes is the number of 1D modes of the expansion \n");
fprintf(stderr," where: - OptimisationLevel is one of the following:\n");
fprintf(stderr," 0: Use elemental sum-factorisation evaluation \n");
fprintf(stderr," 2: Use elemental matrix evaluation using blockmatrices \n");
fprintf(stderr," 3: Use global matrix evaluation \n");
fprintf(stderr," 4: Use optimal evaluation (this option requires optimisation-files being set-up) \n");
exit(1);
}
boost::filesystem::path basePath(BASE_PATH);
int Type = atoi(argv[1]);
int MeshSize = atoi(argv[2]);
int NumModes = atoi(argv[3]);
int optLevel = atoi(argv[4]);
//----------------------------------------------
// Retrieve the necessary input files
stringstream MeshFileName;
stringstream MeshFileDirectory;
stringstream BCfileName;
stringstream ExpansionsFileName;
stringstream GlobOptFileName;
switch(Type)
{
case 1:
{
MeshFileDirectory << "RegularQuadMeshes";
MeshFileName << "UnitSquare_RegularQuadMesh_h_1_" << MeshSize << ".xml";
}
break;
case 2:
{
MeshFileDirectory << "DeformedQuadMeshes";
MeshFileName << "UnitSquare_DeformedQuadMesh_h_1_" << MeshSize << ".xml";
}
break;
case 3:
{
MeshFileDirectory << "RegularTriMeshes";
MeshFileName << "UnitSquare_RegularTriMesh_h_1_" << MeshSize << ".xml";
}
break;
default:
{
cerr << "Type should be equal to one of the following values: "<< endl;
cerr << " 1: Regular Quads" << endl;
cerr << " 2: Deformed Quads" << endl;
cerr << " 3: Regular Tris" << endl;
exit(1);
}
}
BCfileName << "UnitSquare_DirichletBoundaryConditions.xml";
ExpansionsFileName << "NektarExpansionsNummodes" << NumModes << ".xml";
switch(optLevel)
{
case 0:
{
GlobOptFileName << "NoGlobalMat.xml";
}
break;
case 2:
{
GlobOptFileName << "DoBlockMat.xml";
}
break;
case 3:
{
GlobOptFileName << "DoGlobalMat.xml";
}
break;
case 4:
{
ASSERTL0(false,"Optimisation level not set up");
}
break;
default:
{
ASSERTL0(false,"Unrecognised optimisation level");
}
}
boost::filesystem::path MeshFilePath = basePath /
boost::filesystem::path("InputFiles") /
boost::filesystem::path("Geometry") /
boost::filesystem::path(MeshFileDirectory.str()) /
boost::filesystem::path(MeshFileName.str());
vFilenames.push_back(PortablePath(MeshFilePath));
boost::filesystem::path BCfilePath = basePath /
boost::filesystem::path("InputFiles") /
boost::filesystem::path("Conditions") /
boost::filesystem::path(BCfileName.str());
vFilenames.push_back(PortablePath(BCfilePath));
boost::filesystem::path ExpansionsFilePath = basePath /
boost::filesystem::path("InputFiles") /
boost::filesystem::path("Expansions") /
boost::filesystem::path(ExpansionsFileName.str());
vFilenames.push_back(PortablePath(ExpansionsFilePath));
boost::filesystem::path GlobOptFilePath = basePath /
boost::filesystem::path("InputFiles") /
boost::filesystem::path("Optimisation") /
boost::filesystem::path(GlobOptFileName.str());
vFilenames.push_back(PortablePath(GlobOptFilePath));
//----------------------------------------------
LibUtilities::SessionReaderSharedPtr vSession
= LibUtilities::SessionReader::CreateInstance(argc, argv, vFilenames);
//----------------------------------------------
// Read in mesh from input file
SpatialDomains::MeshGraphSharedPtr graph2D = MemoryManager<SpatialDomains::MeshGraph2D>::AllocateSharedPtr(vSession);
//----------------------------------------------
//----------------------------------------------
// Print summary of solution details
lambda = vSession->GetParameter("Lambda");
//----------------------------------------------
//----------------------------------------------
// Define Expansion
Exp = MemoryManager<MultiRegions::ContField2D>::
AllocateSharedPtr(vSession,graph2D,vSession->GetVariable(0));
//----------------------------------------------
int NumElements = Exp->GetExpSize();
//----------------------------------------------
// Set up coordinates of mesh for Forcing function evaluation
coordim = Exp->GetCoordim(0);
nq = Exp->GetTotPoints();
xc0 = Array<OneD,NekDouble>(nq,0.0);
xc1 = Array<OneD,NekDouble>(nq,0.0);
xc2 = Array<OneD,NekDouble>(nq,0.0);
switch(coordim)
{
case 1:
Exp->GetCoords(xc0);
break;
case 2:
Exp->GetCoords(xc0,xc1);
break;
case 3:
Exp->GetCoords(xc0,xc1,xc2);
break;
}
//----------------------------------------------
//----------------------------------------------
// Define forcing function for first variable defined in file
fce = Array<OneD,NekDouble>(nq);
LibUtilities::EquationSharedPtr ffunc = vSession->GetFunction("Forcing",0);
for(i = 0; i < nq; ++i)
{
fce[i] = ffunc->Evaluate(xc0[i],xc1[i],xc2[i]);
}
//----------------------------------------------
//----------------------------------------------
// Setup expansion containing the forcing function
Fce = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(*Exp);
Fce->SetPhys(fce);
//----------------------------------------------
//----------------------------------------------
// Helmholtz solution taking physical forcing
FlagList flags;
flags.set(eUseGlobal, true);
StdRegions::ConstFactorMap factors;
factors[StdRegions::eFactorLambda] = lambda;
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(),flags,factors);
//----------------------------------------------
//----------------------------------------------
// Backward Transform Solution to get solved values at
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(),
MultiRegions::eGlobal);
//----------------------------------------------
//----------------------------------------------
// See if there is an exact solution, if so
// evaluate and plot errors
LibUtilities::EquationSharedPtr ex_sol = vSession->GetFunction("ExactSolution",0);
//----------------------------------------------
// evaluate exact solution
sol = Array<OneD,NekDouble>(nq);
for(i = 0; i < nq; ++i)
{
sol[i] = ex_sol->Evaluate(xc0[i],xc1[i],xc2[i]);
}
//----------------------------------------------
//--------------------------------------------
// Calculate L_inf error
Sol = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(*Exp);
Sol->SetPhys(sol);
Sol->SetPhysState(true);
NekDouble L2Error = Exp->L2 (Sol->GetPhys());
NekDouble LinfError = Exp->Linf(Sol->GetPhys());
//--------------------------------------------
// alternative error calculation
const LibUtilities::PointsKey PkeyT1(30,LibUtilities::eGaussLobattoLegendre);
const LibUtilities::PointsKey PkeyT2(30,LibUtilities::eGaussRadauMAlpha1Beta0);
const LibUtilities::PointsKey PkeyQ1(30,LibUtilities::eGaussLobattoLegendre);
const LibUtilities::PointsKey PkeyQ2(30,LibUtilities::eGaussLobattoLegendre);
const LibUtilities::BasisKey BkeyT1(LibUtilities::eModified_A,NumModes,PkeyT1);
const LibUtilities::BasisKey BkeyT2(LibUtilities::eModified_B,NumModes,PkeyT2);
const LibUtilities::BasisKey BkeyQ1(LibUtilities::eModified_A,NumModes,PkeyQ1);
const LibUtilities::BasisKey BkeyQ2(LibUtilities::eModified_A,NumModes,PkeyQ2);
MultiRegions::ExpList2DSharedPtr ErrorExp =
MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,BkeyT1,BkeyT2,BkeyQ1,BkeyQ2,graph2D);
int ErrorCoordim = ErrorExp->GetCoordim(0);
int ErrorNq = ErrorExp->GetTotPoints();
Array<OneD,NekDouble> ErrorXc0(ErrorNq,0.0);
Array<OneD,NekDouble> ErrorXc1(ErrorNq,0.0);
Array<OneD,NekDouble> ErrorXc2(ErrorNq,0.0);
switch(ErrorCoordim)
{
case 1:
ErrorExp->GetCoords(ErrorXc0);
break;
case 2:
ErrorExp->GetCoords(ErrorXc0,ErrorXc1);
break;
case 3:
ErrorExp->GetCoords(ErrorXc0,ErrorXc1,ErrorXc2);
break;
}
// evaluate exact solution
Array<OneD,NekDouble> ErrorSol(ErrorNq);
for(i = 0; i < ErrorNq; ++i)
{
ErrorSol[i] = ex_sol->Evaluate(ErrorXc0[i],ErrorXc1[i],ErrorXc2[i]);
}
// calcualte spectral/hp approximation on the quad points of this new
// expansion basis
Exp->GlobalToLocal(Exp->GetCoeffs(),ErrorExp->UpdateCoeffs());
ErrorExp->BwdTrans_IterPerExp(ErrorExp->GetCoeffs(),ErrorExp->UpdatePhys());
NekDouble L2ErrorBis = ErrorExp->L2 (ErrorSol);
NekDouble LinfErrorBis = ErrorExp->Linf(ErrorSol);
//--------------------------------------------
MultiRegions::ContField2DSharedPtr Exp,Fce,Sol;
int i, nq, coordim;
Array<OneD,NekDouble> fce,sol;
Array<OneD,NekDouble> xc0,xc1,xc2;
NekDouble lambda;
vector<string> vFilenames;
if(argc != 5)
{
fprintf(stderr,"Usage: TimingCGHelmSolve2D Type MeshSize NumModes OptimisationLevel\n");
fprintf(stderr," where: - Type is one of the following:\n");
fprintf(stderr," 1: Regular Quadrilaterals \n");
fprintf(stderr," 2: Deformed Quadrilaterals (may not be supported) \n");
fprintf(stderr," 3: Regular Triangles \n");
fprintf(stderr," where: - MeshSize is 1/h \n");
fprintf(stderr," where: - NumModes is the number of 1D modes of the expansion \n");
fprintf(stderr," where: - OptimisationLevel is one of the following:\n");
fprintf(stderr," 0: Use elemental sum-factorisation evaluation \n");
fprintf(stderr," 2: Use elemental matrix evaluation using blockmatrices \n");
fprintf(stderr," 3: Use global matrix evaluation \n");
fprintf(stderr," 4: Use optimal evaluation (this option requires optimisation-files being set-up) \n");
exit(1);
}
boost::filesystem::path basePath(BASE_PATH);
int Type = atoi(argv[1]);
std::string TypeStr;
int MeshSize = atoi(argv[2]);
int NumModes = atoi(argv[3]);
int optLevel = atoi(argv[4]);
std::string optLevelStr;
//----------------------------------------------
// Retrieve the necessary input files
stringstream MeshFileName;
stringstream MeshFileDirectory;
stringstream BCfileName;
stringstream ExpansionsFileName;
stringstream GlobOptFileName;
switch(Type)
{
case 1:
{
MeshFileDirectory << "RegularQuadMeshes";
MeshFileName << "UnitSquare_RegularQuadMesh_h_1_" << MeshSize << ".xml";
TypeStr = "RQuad";
}
break;
case 2:
{
MeshFileDirectory << "DeformedQuadMeshes";
MeshFileName << "UnitSquare_DeformedQuadMesh_h_1_" << MeshSize << ".xml";
TypeStr = "DQuad";
}
break;
case 3:
{
MeshFileDirectory << "RegularTriMeshes";
MeshFileName << "UnitSquare_RegularTriMesh_h_1_" << MeshSize << ".xml";
TypeStr = "RTri";
}
break;
default:
{
cerr << "Type should be equal to one of the following values: "<< endl;
cerr << " 1: Regular Quads" << endl;
cerr << " 2: Deformed Quads" << endl;
cerr << " 3: Regular Tris" << endl;
exit(1);
}
}
BCfileName << "UnitSquare_DirichletBoundaryConditions.xml";
ExpansionsFileName << "NektarExpansionsNummodes" << NumModes << ".xml";
switch(optLevel)
{
case 0:
{
GlobOptFileName << "NoGlobalMat.xml";
optLevelStr = "SumFac";
}
break;
case 2:
{
GlobOptFileName << "DoBlockMat.xml";
optLevelStr = "BlkMat";
}
break;
case 3:
{
GlobOptFileName << "DoGlobalMat.xml";
optLevelStr = "GlbMat";
}
break;
case 4:
{
ASSERTL0(false,"Optimisation level not set up");
}
break;
default:
{
ASSERTL0(false,"Unrecognised optimisation level");
}
}
boost::filesystem::path MeshFilePath = basePath /
boost::filesystem::path("InputFiles") /
boost::filesystem::path("Geometry") /
boost::filesystem::path(MeshFileDirectory.str()) /
boost::filesystem::path(MeshFileName.str());
vFilenames.push_back(PortablePath(MeshFilePath));
boost::filesystem::path BCfilePath = basePath /
boost::filesystem::path("InputFiles") /
boost::filesystem::path("Conditions") /
boost::filesystem::path(BCfileName.str());
vFilenames.push_back(PortablePath(BCfilePath));
boost::filesystem::path ExpansionsFilePath = basePath /
boost::filesystem::path("InputFiles") /
boost::filesystem::path("Expansions") /
boost::filesystem::path(ExpansionsFileName.str());
vFilenames.push_back(PortablePath(ExpansionsFilePath));
boost::filesystem::path GlobOptFilePath = basePath /
boost::filesystem::path("InputFiles") /
boost::filesystem::path("Optimisation") /
boost::filesystem::path(GlobOptFileName.str());
vFilenames.push_back(PortablePath(GlobOptFilePath));
//----------------------------------------------
LibUtilities::SessionReaderSharedPtr vSession
= LibUtilities::SessionReader::CreateInstance(argc, argv, vFilenames);
//----------------------------------------------
// Read in mesh from input file
SpatialDomains::MeshGraphSharedPtr graph2D = MemoryManager<SpatialDomains::MeshGraph2D>::AllocateSharedPtr(vSession);
//----------------------------------------------
//----------------------------------------------
// Print summary of solution details
lambda = vSession->GetParameter("Lambda");
//----------------------------------------------
//----------------------------------------------
// Define Expansion
Exp = MemoryManager<MultiRegions::ContField2D>::
AllocateSharedPtr(vSession,graph2D,vSession->GetVariable(0));
//----------------------------------------------
int NumElements = Exp->GetExpSize();
//----------------------------------------------
// Set up coordinates of mesh for Forcing function evaluation
coordim = Exp->GetCoordim(0);
nq = Exp->GetTotPoints();
xc0 = Array<OneD,NekDouble>(nq,0.0);
xc1 = Array<OneD,NekDouble>(nq,0.0);
xc2 = Array<OneD,NekDouble>(nq,0.0);
switch(coordim)
{
case 1:
Exp->GetCoords(xc0);
break;
case 2:
Exp->GetCoords(xc0,xc1);
break;
case 3:
Exp->GetCoords(xc0,xc1,xc2);
break;
}
//----------------------------------------------
//----------------------------------------------
// Define forcing function for first variable defined in file
fce = Array<OneD,NekDouble>(nq);
LibUtilities::EquationSharedPtr ffunc = vSession->GetFunction("Forcing",0);
ffunc->Evaluate(xc0,xc1,xc2,fce);
//----------------------------------------------
//----------------------------------------------
// Setup expansion containing the forcing function
Fce = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(*Exp);
Fce->SetPhys(fce);
//----------------------------------------------
//----------------------------------------------
// Helmholtz solution taking physical forcing
FlagList flags;
flags.set(eUseGlobal, true);
StdRegions::ConstFactorMap factors;
factors[StdRegions::eFactorLambda] = lambda;
Exp->HelmSolve(Fce->GetPhys(), Exp->UpdateCoeffs(),flags,factors);
//----------------------------------------------
//----------------------------------------------
// Backward Transform Solution to get solved values at
Exp->BwdTrans(Exp->GetCoeffs(), Exp->UpdatePhys(),
MultiRegions::eGlobal);
//----------------------------------------------
//----------------------------------------------
// See if there is an exact solution, if so
// evaluate and plot errors
LibUtilities::EquationSharedPtr ex_sol = vSession->GetFunction("ExactSolution",0);
//----------------------------------------------
// evaluate exact solution
sol = Array<OneD,NekDouble>(nq);
ex_sol->Evaluate(xc0,xc1,xc2,sol);
//----------------------------------------------
//--------------------------------------------
// Calculate L_inf error
Sol = MemoryManager<MultiRegions::ContField2D>::AllocateSharedPtr(*Exp);
NekDouble L2Error = Exp->L2 (sol);
NekDouble LinfError = Exp->Linf(sol);
//--------------------------------------------
// alternative error calculation
const LibUtilities::PointsKey PkeyT1(30,LibUtilities::eGaussLobattoLegendre);
const LibUtilities::PointsKey PkeyT2(30,LibUtilities::eGaussRadauMAlpha1Beta0);
const LibUtilities::PointsKey PkeyQ1(30,LibUtilities::eGaussLobattoLegendre);
const LibUtilities::PointsKey PkeyQ2(30,LibUtilities::eGaussLobattoLegendre);
const LibUtilities::BasisKey BkeyT1(LibUtilities::eModified_A,NumModes,PkeyT1);
const LibUtilities::BasisKey BkeyT2(LibUtilities::eModified_B,NumModes,PkeyT2);
const LibUtilities::BasisKey BkeyQ1(LibUtilities::eModified_A,NumModes,PkeyQ1);
const LibUtilities::BasisKey BkeyQ2(LibUtilities::eModified_A,NumModes,PkeyQ2);
MultiRegions::ExpList2DSharedPtr ErrorExp =
MemoryManager<MultiRegions::ExpList2D>::AllocateSharedPtr(vSession,BkeyT1,BkeyT2,BkeyQ1,BkeyQ2,graph2D);
int ErrorCoordim = ErrorExp->GetCoordim(0);
int ErrorNq = ErrorExp->GetTotPoints();
Array<OneD,NekDouble> ErrorXc0(ErrorNq,0.0);
Array<OneD,NekDouble> ErrorXc1(ErrorNq,0.0);
Array<OneD,NekDouble> ErrorXc2(ErrorNq,0.0);
switch(ErrorCoordim)