Nektar++
Aliasing.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstdlib>
3 
5 #include <SolverUtils/Driver.h>
6 
9 
10 using namespace std;
11 using namespace Nektar;
12 using namespace Nektar::SolverUtils;
13 
14 int main(int argc, char *argv[])
15 {
16  if (argc != 2)
17  {
18  fprintf(stderr, "Usage: ./Aliasing file.xml \n");
19  fprintf(stderr, "\t Method will read intiial conditions section of "
20  ".xml file for input \n");
21  exit(1);
22  }
23 
26  string vDriverModule;
27  DriverSharedPtr drv;
28  try
29  {
30  // Create session reader.
31  session = LibUtilities::SessionReader::CreateInstance(argc, argv);
32 
33  // Create MeshGraph.
34  graph = SpatialDomains::MeshGraph::Read(session);
35 
36  // Create driver
37  session->LoadSolverInfo("Driver", vDriverModule, "Standard");
38  drv = GetDriverFactory().CreateInstance(vDriverModule, session, graph);
39 
40  EquationSystemSharedPtr EqSys = drv->GetEqu()[0];
41  IncNavierStokesSharedPtr IncNav = EqSys->as<IncNavierStokes>();
42 
43  IncNav->SetInitialConditions(0.0, false);
45  IncNav->UpdateFields();
46 
47  int i;
48  int nConvectiveFields = IncNav->GetNConvectiveFields();
49  int nphys = fields[0]->GetTotPoints();
50  Array<OneD, Array<OneD, NekDouble>> VelFields(nConvectiveFields);
51  Array<OneD, Array<OneD, NekDouble>> NonLinear(nConvectiveFields);
52  Array<OneD, Array<OneD, NekDouble>> NonLinearDealiased(
53  nConvectiveFields);
54 
55  for (i = 0; i < nConvectiveFields; ++i)
56  {
57  VelFields[i] = fields[i]->UpdatePhys();
58  NonLinear[i] = Array<OneD, NekDouble>(nphys);
59  NonLinearDealiased[i] = Array<OneD, NekDouble>(nphys);
60  }
61 
62  std::shared_ptr<NavierStokesAdvection> A =
63  std::dynamic_pointer_cast<NavierStokesAdvection>(
64  IncNav->GetAdvObject());
65 
66  if (!A)
67  {
68  cout << "Must use non-linear Navier-Stokes advection" << endl;
69  exit(-1);
70  }
71 
72  // calculate non-linear terms without dealiasing
73  A->SetSpecHPDealiasing(false);
74  A->Advect(nConvectiveFields, fields, VelFields, VelFields, NonLinear,
75  0.0);
76 
77  // calculate non-linear terms with dealiasing
78  A->SetSpecHPDealiasing(true);
79  A->Advect(nConvectiveFields, fields, VelFields, VelFields,
80  NonLinearDealiased, 0.0);
81 
82  // Evaulate Difference and put into fields;
83  for (i = 0; i < nConvectiveFields; ++i)
84  {
85  Vmath::Vsub(nphys, NonLinearDealiased[i], 1, NonLinear[i], 1,
86  NonLinear[i], 1);
87  fields[i]->FwdTransLocalElmt(NonLinear[i],
88  fields[i]->UpdateCoeffs());
89  // Need to reset varibale name for output
90  string name = "NL_Aliasing_" + session->GetVariable(i);
91  session->SetVariable(i, name.c_str());
92  }
93 
94  // Reset session name for output file
95  std::string outname = IncNav->GetSessionName();
96 
97  outname += "_NonLinear_Aliasing";
98  IncNav->ResetSessionName(outname);
99  IncNav->Output();
100  }
101  catch (const std::runtime_error &)
102  {
103  return 1;
104  }
105  catch (const std::string &eStr)
106  {
107  cout << "Error: " << eStr << endl;
108  }
109 
110  return 0;
111 }
int main(int argc, char *argv[])
Definition: Aliasing.cpp:14
This class is the base class for Navier Stokes problems.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
SOLVER_UTILS_EXPORT void SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Initialise the data in the dependent fields.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:51
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
std::shared_ptr< EquationSystem > EquationSystemSharedPtr
A shared pointer to an EquationSystem object.
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< IncNavierStokes > IncNavierStokesSharedPtr
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.
Definition: Vmath.cpp:419