Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
DiffusionSolver.cpp File Reference
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <SpatialDomains/MeshGraph.h>
#include <MultiRegions/ContField2D.h>
Include dependency graph for DiffusionSolver.cpp:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 43 of file DiffusionSolver.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), Nektar::LibUtilities::SessionReader::CreateInstance(), Nektar::StdRegions::eFactorLambda, Nektar::NullFlagList, Nektar::SpatialDomains::MeshGraph::Read(), Vmath::Smul(), and Vmath::Zero().

{
try
{
// Create session reader.
session = LibUtilities::SessionReader::CreateInstance(argc, argv);
// Create Field I/O object.
AllocateSharedPtr(session->GetComm());
// Get some information about the session
string sessionName = session->GetSessionName();
string outFile = sessionName + ".fld";
unsigned int nSteps = session->GetParameter("NumSteps");
NekDouble delta_t = session->GetParameter("TimeStep");
NekDouble epsilon = session->GetParameter("epsilon" );
// Read the geometry and the expansion information
graph = SpatialDomains::MeshGraph::Read(session);
// Create field
::AllocateSharedPtr(session, graph, session->GetVariable(0));
// Get coordinates of physical points
unsigned int nq = field->GetNpoints();
Array<OneD,NekDouble> x0(nq), x1(nq), x2(nq);
field->GetCoords(x0,x1,x2);
// Evaluate initial condition at these points
icond = session->GetFunction("InitialConditions", "u");
icond->Evaluate(x0, x1, x2, 0.0, field->UpdatePhys());
// Compute lambda in the Helmholtz problem
factors[StdRegions::eFactorLambda] = 1.0/delta_t/epsilon;
// Zero field coefficients for initial guess for linear solver.
Vmath::Zero(field->GetNcoeffs(), field->UpdateCoeffs(), 1);
// Time integrate using backward Euler
for (unsigned int n = 0; n < nSteps; ++n)
{
Vmath::Smul(nq, -1.0/delta_t/epsilon, field->GetPhys(), 1,
field->UpdatePhys(), 1);
field->HelmSolve(field->GetPhys(), field->UpdateCoeffs(),
NullFlagList, factors);
field->BwdTrans(field->GetCoeffs(), field->UpdatePhys());
}
// Write solution to file
std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
= field->GetFieldDefinitions();
std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
for(int i = 0; i < FieldDef.size(); ++i)
{
FieldDef[i]->m_fields.push_back("u");
field->AppendFieldData(FieldDef[i], FieldData[i]);
}
fld->Write(outFile, FieldDef, FieldData);
// Check for exact solution
ex_sol = session->GetFunction("ExactSolution",0);
if(ex_sol)
{
// Allocate storage
Array<OneD, NekDouble> exact(nq);
//----------------------------------------------
// evaluate exact solution
ex_sol->Evaluate(x0, x1, x2, (nSteps)*delta_t, exact);
//--------------------------------------------
// Calculate errors
cout << "L inf error: "
<< field->Linf(field->GetPhys(), exact) << endl;
cout << "L 2 error: "
<< field->L2(field->GetPhys(), exact) << endl;
cout << "H 1 error: "
<< field->H1(field->GetPhys(), exact) << endl;
//--------------------------------------------
}
// Finalise session
session->Finalise();
}
catch (const std::runtime_error& e)
{
return 1;
}
catch (const std::string& eStr)
{
cout << "Error: " << eStr << endl;
}
return 0;
}