Nektar++
DiffusionSolver.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ADRSolver.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Advection Diffusion Reaction framework solver
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
40 
41 using namespace Nektar;
42 
43 int main(int argc, char *argv[])
44 {
49  LibUtilities::EquationSharedPtr icond, ex_sol;
51 
52  try
53  {
54  // Create session reader.
55  session = LibUtilities::SessionReader::CreateInstance(argc, argv);
56 
57  // Create Field I/O object.
59  AllocateSharedPtr(session->GetComm());
60 
61  // Get some information about the session
62  string sessionName = session->GetSessionName();
63  string outFile = sessionName + ".fld";
64  unsigned int nSteps = session->GetParameter("NumSteps");
65  NekDouble delta_t = session->GetParameter("TimeStep");
66  NekDouble epsilon = session->GetParameter("epsilon" );
67 
68  // Read the geometry and the expansion information
69  graph = SpatialDomains::MeshGraph::Read(session);
70 
71  // Create field
73  ::AllocateSharedPtr(session, graph, session->GetVariable(0));
74 
75  // Get coordinates of physical points
76  unsigned int nq = field->GetNpoints();
77  Array<OneD,NekDouble> x0(nq), x1(nq), x2(nq);
78  field->GetCoords(x0,x1,x2);
79 
80  // Evaluate initial condition at these points
81  icond = session->GetFunction("InitialConditions", "u");
82  icond->Evaluate(x0, x1, x2, 0.0, field->UpdatePhys());
83 
84  // Compute lambda in the Helmholtz problem
85  factors[StdRegions::eFactorLambda] = 1.0/delta_t/epsilon;
86 
87  // Zero field coefficients for initial guess for linear solver.
88  Vmath::Zero(field->GetNcoeffs(), field->UpdateCoeffs(), 1);
89 
90  // Time integrate using backward Euler
91  for (unsigned int n = 0; n < nSteps; ++n)
92  {
93  Vmath::Smul(nq, -1.0/delta_t/epsilon, field->GetPhys(), 1,
94  field->UpdatePhys(), 1);
95 
96  field->HelmSolve(field->GetPhys(), field->UpdateCoeffs(),
97  NullFlagList, factors);
98 
99  field->BwdTrans(field->GetCoeffs(), field->UpdatePhys());
100  }
101 
102  // Write solution to file
103  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
104  = field->GetFieldDefinitions();
105  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
106  for(int i = 0; i < FieldDef.size(); ++i)
107  {
108  FieldDef[i]->m_fields.push_back("u");
109  field->AppendFieldData(FieldDef[i], FieldData[i]);
110  }
111  fld->Write(outFile, FieldDef, FieldData);
112 
113  // Check for exact solution
114  ex_sol = session->GetFunction("ExactSolution",0);
115  if(ex_sol)
116  {
117  // Allocate storage
118  Array<OneD, NekDouble> exact(nq);
119 
120  //----------------------------------------------
121  // evaluate exact solution
122  ex_sol->Evaluate(x0, x1, x2, (nSteps)*delta_t, exact);
123 
124  //--------------------------------------------
125  // Calculate errors
126  cout << "L inf error: "
127  << field->Linf(field->GetPhys(), exact) << endl;
128  cout << "L 2 error: "
129  << field->L2(field->GetPhys(), exact) << endl;
130  cout << "H 1 error: "
131  << field->H1(field->GetPhys(), exact) << endl;
132  //--------------------------------------------
133  }
134 
135  // Finalise session
136  session->Finalise();
137  }
138  catch (const std::runtime_error& e)
139  {
140  return 1;
141  }
142  catch (const std::string& eStr)
143  {
144  cout << "Error: " << eStr << endl;
145  }
146 
147  return 0;
148 }
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:119
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:291
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:236
double NekDouble
int main(int argc, char *argv[])
boost::shared_ptr< Equation > EquationSharedPtr
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
static FlagList NullFlagList
An empty flag list.