Commit de03bc73 authored by Chris Cantwell's avatar Chris Cantwell

Refactored compressible flow solver main function.

Generalised error calculation and output.
Linf now computes max value in solution, rather than zero, to be consistent
with L2 equivalent.
parent a68fd172
...@@ -33,6 +33,8 @@ ...@@ -33,6 +33,8 @@
// //
/////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <SolverUtils/DriverStandard.h> #include <SolverUtils/DriverStandard.h>
namespace Nektar namespace Nektar
...@@ -71,23 +73,49 @@ namespace Nektar ...@@ -71,23 +73,49 @@ namespace Nektar
void DriverStandard::v_Execute(ostream &out) void DriverStandard::v_Execute(ostream &out)
{ {
time_t starttime, endtime;
NekDouble CPUtime;
m_equ[0]->PrintSummary(out); m_equ[0]->PrintSummary(out);
time(&starttime);
m_equ[0]->DoInitialise(); m_equ[0]->DoInitialise();
m_equ[0]->DoSolve(); m_equ[0]->DoSolve();
time(&endtime);
m_equ[0]->Output(); m_equ[0]->Output();
if (m_comm->GetRank() == 0)
{
CPUtime = difftime(endtime, starttime);
cout << "-------------------------------------------" << endl;
cout << "Total Computation Time = " << CPUtime << "s" << endl;
cout << "-------------------------------------------" << endl;
}
// Evaluate and output computation time and solution accuracy. // Evaluate and output computation time and solution accuracy.
// The specific format of the error output is essential for the // The specific format of the error output is essential for the
// regression tests to work. // regression tests to work.
// Evaluate L2 Error // Evaluate L2 Error
for(int i = 0; i < m_equ[0]->GetNvariables(); ++i) for(int i = 0; i < m_equ[0]->GetNvariables(); ++i)
{ {
NekDouble vL2Error = m_equ[0]->L2Error(i,false); Array<OneD, NekDouble> exactsoln(m_equ[0]->GetTotPoints(), 0.0);
NekDouble vLinfError = m_equ[0]->LinfError(i);
// Evaluate "ExactSolution" function, or zero array
m_equ[0]->EvaluateExactSolution(i, exactsoln,
m_equ[0]->GetFinalTime());
NekDouble vL2Error = m_equ[0]->L2Error (i, exactsoln);
NekDouble vLinfError = m_equ[0]->LinfError(i, exactsoln);
if (m_comm->GetRank() == 0) if (m_comm->GetRank() == 0)
{ {
out << "L 2 error (variable " << m_equ[0]->GetVariable(i) << ") : " << vL2Error << endl; out << "L 2 error (variable " << m_equ[0]->GetVariable(i)
out << "L inf error (variable " << m_equ[0]->GetVariable(i) << ") : " << vLinfError << endl; << ") : " << vL2Error << endl;
out << "L inf error (variable " << m_equ[0]->GetVariable(i)
<< ") : " << vLinfError << endl;
} }
} }
} }
......
...@@ -33,7 +33,6 @@ ...@@ -33,7 +33,6 @@
// //
/////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////
#include <SolverUtils/EquationSystem.h>
#include <SolverUtils/DriverSteadyState.h> #include <SolverUtils/DriverSteadyState.h>
...@@ -76,26 +75,17 @@ namespace Nektar ...@@ -76,26 +75,17 @@ namespace Nektar
//With a loop over "DoSolve", this Driver implements the "encaplulated" Selective Frequency Damping method //With a loop over "DoSolve", this Driver implements the "encaplulated" Selective Frequency Damping method
//to find the steady state of a flow above the critical Reynolds number. //to find the steady state of a flow above the critical Reynolds number.
//Read the graph to define the dimension of the problem
SpatialDomains::MeshGraphSharedPtr graph;
graph = SpatialDomains::MeshGraph::Read(m_session);
m_equ[0]->PrintSummary(out); m_equ[0]->PrintSummary(out);
m_equ[0]->DoInitialise(); m_equ[0]->DoInitialise();
// - SFD Routine - // - SFD Routine -
// Compressible case
if (m_session->GetSolverInfo("EqType") == "EulerCFE" || m_session->GetSolverInfo("EqType") == "NavierStokesCFE") //Compressible case NumVar_SFD = m_equ[0]->UpdateFields()[0]->GetCoordim(0);
if (m_session->GetSolverInfo("EqType") == "EulerCFE" ||
m_session->GetSolverInfo("EqType") == "NavierStokesCFE")
{ {
cout << "Compressible case!!!" << endl; NumVar_SFD += 2; //Number of variables for the compressible equations
NumVar_SFD = graph->GetSpaceDimension() + 2; //Number of variables for the compressible equations
} }
else //Incompressible case
{
cout << "InCompressible case!!!" << endl;
NumVar_SFD = graph->GetSpaceDimension(); //Number of velocity components for the incompressible equations
}
Array<OneD, Array<OneD, NekDouble> > q0(NumVar_SFD); Array<OneD, Array<OneD, NekDouble> > q0(NumVar_SFD);
Array<OneD, Array<OneD, NekDouble> > q1(NumVar_SFD); Array<OneD, Array<OneD, NekDouble> > q1(NumVar_SFD);
......
...@@ -963,7 +963,7 @@ namespace Nektar ...@@ -963,7 +963,7 @@ namespace Nektar
} }
else else
{ {
Linferror = 0.0; Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys());
} }
} }
else else
...@@ -1109,13 +1109,14 @@ namespace Nektar ...@@ -1109,13 +1109,14 @@ namespace Nektar
Array<OneD, NekDouble> &outfield, Array<OneD, NekDouble> &outfield,
const NekDouble time) const NekDouble time)
{ {
ASSERTL0 (m_session->DefinesFunction("ExactSolution"),
"No ExactSolution provided in session file.");
ASSERTL0 (outfield.num_elements() == m_fields[field]->GetNpoints(), ASSERTL0 (outfield.num_elements() == m_fields[field]->GetNpoints(),
"ExactSolution array size mismatch."); "ExactSolution array size mismatch.");
Vmath::Zero(outfield.num_elements(), outfield, 1);
EvaluateFunction(m_session->GetVariable(field), outfield, if (m_session->DefinesFunction("ExactSolution"))
"ExactSolution", time); {
EvaluateFunction(m_session->GetVariable(field), outfield,
"ExactSolution", time);
}
} }
......
...@@ -42,84 +42,34 @@ using namespace Nektar::SolverUtils; ...@@ -42,84 +42,34 @@ using namespace Nektar::SolverUtils;
int main(int argc, char *argv[]) int main(int argc, char *argv[])
{ {
// Create session reader
LibUtilities::SessionReaderSharedPtr session; LibUtilities::SessionReaderSharedPtr session;
session = LibUtilities::SessionReader::CreateInstance(argc, argv); string vDriverModule;
DriverSharedPtr drv;
// === Added for using SteadyState driver for compressible solver
if (session->GetSolverInfo("Driver") == "SteadyState")
{
string vDriverModule;
DriverSharedPtr drv;
session->LoadSolverInfo("Driver", vDriverModule, "SteadyState");
drv = GetDriverFactory().CreateInstance(vDriverModule, session);
drv->Execute();
}
// === end ===
time_t starttime, endtime;
NekDouble CPUtime;
EquationSystemSharedPtr equ;
// Record start time
time(&starttime);
// Create instance of module to solve the equation specified in the session
try try
{ {
equ = GetEquationSystemFactory().CreateInstance( // Create session reader.
session->GetSolverInfo("EQTYPE"), session); session = LibUtilities::SessionReader::CreateInstance(argc, argv);
}
catch (int e)
{
ASSERTL0(e == -1, "No such solver class defined.");
}
// Print a summary of solver/problem parameters
equ->PrintSummary(cout);
// Initialise the problem
equ->DoInitialise();
// Solve the problem // Create driver
equ->DoSolve(); session->LoadSolverInfo("Driver", vDriverModule, "Standard");
drv = GetDriverFactory().CreateInstance(vDriverModule, session);
// Record end time
time(&endtime);
// Compute the computational time in hours
CPUtime = difftime(endtime, starttime);
// Write output to .fld file // Execute driver
equ->Output(); drv->Execute();
// Evaluate and output computation time and solution accuracy // Finalise communications
if (session->GetComm()->GetRank() == 0) session->Finalise();
}
catch (const std::runtime_error&)
{ {
cout << "-------------------------------------------" << endl; return 1;
cout << "Total Computation Time = " << CPUtime << "s" << endl;
cout << "-------------------------------------------" << endl;
} }
catch (const std::string& eStr)
for(int i = 0; i < equ->GetNvariables(); ++i)
{ {
// Get exact solution cout << "Error: " << eStr << endl;
Array<OneD, NekDouble> exactsoln(equ->GetTotPoints(), 0.0);
equ->EvaluateExactSolution(i, exactsoln, equ->GetFinalTime());
NekDouble l2 = equ->L2Error (i, exactsoln);
NekDouble li = equ->LinfError(i, exactsoln);
if (session->GetComm()->GetRank() == 0)
{
cout << "L 2 error (variable " << equ->GetVariable(i) << "): "
<< l2 << endl;
cout << "L inf error (variable " << equ->GetVariable(i) << "): "
<< li << endl;
}
} }
session->Finalise(); return 0;
} }
...@@ -15,9 +15,9 @@ ...@@ -15,9 +15,9 @@
<value variable="p" tolerance="1e-12">0.0123609</value> <value variable="p" tolerance="1e-12">0.0123609</value>
</metric> </metric>
<metric type="Linf" id="2"> <metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">0</value> <value variable="u" tolerance="1e-12">0.0272836</value>
<value variable="v" tolerance="1e-12">0</value> <value variable="v" tolerance="1e-12">0.0122607</value>
<value variable="p" tolerance="1e-12">0</value> <value variable="p" tolerance="1e-12">0.00217552</value>
</metric> </metric>
</metrics> </metrics>
</test> </test>
......
...@@ -13,9 +13,9 @@ ...@@ -13,9 +13,9 @@
<value variable="p" tolerance="1e-12">0.000232476</value> <value variable="p" tolerance="1e-12">0.000232476</value>
</metric> </metric>
<metric type="Linf" id="2"> <metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">0</value> <value variable="u" tolerance="1e-12">0.000924245</value>
<value variable="v" tolerance="1e-12">0</value> <value variable="v" tolerance="1e-12">0.000576232</value>
<value variable="p" tolerance="1e-12">0</value> <value variable="p" tolerance="1e-12">0.000671216</value>
</metric> </metric>
</metrics> </metrics>
</test> </test>
......
...@@ -13,9 +13,9 @@ ...@@ -13,9 +13,9 @@
<value variable="p" tolerance="1e-12">0.0128673</value> <value variable="p" tolerance="1e-12">0.0128673</value>
</metric> </metric>
<metric type="Linf" id="2"> <metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">0</value> <value variable="u" tolerance="1e-12">1</value>
<value variable="v" tolerance="1e-12">0</value> <value variable="v" tolerance="1e-12">0.000630178</value>
<value variable="p" tolerance="1e-12">0</value> <value variable="p" tolerance="1e-12">0.00746041</value>
</metric> </metric>
</metrics> </metrics>
</test> </test>
......
...@@ -15,9 +15,9 @@ ...@@ -15,9 +15,9 @@
<value variable="p" tolerance="1e-12">0.0157668</value> <value variable="p" tolerance="1e-12">0.0157668</value>
</metric> </metric>
<metric type="Linf" id="2"> <metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">0</value> <value variable="u" tolerance="1e-12">0.0434598</value>
<value variable="v" tolerance="1e-12">0</value> <value variable="v" tolerance="1e-12">0.0467868</value>
<value variable="p" tolerance="1e-12">0</value> <value variable="p" tolerance="1e-12">0.0232832</value>
</metric> </metric>
</metrics> </metrics>
</test> </test>
......
...@@ -50,7 +50,7 @@ namespace Nektar ...@@ -50,7 +50,7 @@ namespace Nektar
// name if it exists: first field is variable name, second field is L2 // name if it exists: first field is variable name, second field is L2
// error. // error.
m_regex = "^L 2 error\\s*(?:\\(variable " m_regex = "^L 2 error\\s*(?:\\(variable "
"(\\w+)\\))?\\s*:\\s*([+-]?\\d.+\\d|-?0|[+-]?nan|[+-]?inf).*"; "(\\w+)\\))?\\s*:\\s*([+-]?\\d.+\\d|-?\\d|[+-]?nan|[+-]?inf).*";
// Find the L2 error to match against. // Find the L2 error to match against.
TiXmlElement *value = metric->FirstChildElement("value"); TiXmlElement *value = metric->FirstChildElement("value");
......
...@@ -50,7 +50,7 @@ namespace Nektar ...@@ -50,7 +50,7 @@ namespace Nektar
// name if it exists: first field is variable name, second field is L2 // name if it exists: first field is variable name, second field is L2
// error. // error.
m_regex = "^L inf\\w* error\\s*(?:\\(variable " m_regex = "^L inf\\w* error\\s*(?:\\(variable "
"(\\w+)\\))?\\s*:\\s*([+-]?\\d.+\\d|-?0|[+-]?nan|[+-]?inf).*"; "(\\w+)\\))?\\s*:\\s*([+-]?\\d.+\\d|-?\\d|[+-]?nan|[+-]?inf).*";
// Find the L2 error to match against. // Find the L2 error to match against.
TiXmlElement *value = metric->FirstChildElement("value"); TiXmlElement *value = metric->FirstChildElement("value");
......
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