Nektar++
DiffusionSolverTimeInt.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File DiffusionTestTI.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: Diffusion solver
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <cstdlib>
37 
43 
44 using namespace Nektar;
45 
46 class Diffusion
47 {
48  public:
49  Diffusion(int argc, char* argv[]);
50  ~Diffusion();
51 
52  void TimeIntegrate();
53 
54  void DoImplicitSolve(
55  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
56  Array<OneD, Array<OneD, NekDouble> >&outarray,
57  const NekDouble time,
58  const NekDouble lambda);
59 
60  private:
63  string sessionName;
66 
71 
72  string scheme;
73  unsigned int nSteps;
77 
78  void WriteSolution();
79  void ExactSolution();
80 };
81 
82 
83 Diffusion::Diffusion(int argc, char* argv[])
84 {
85  // Create session reader.
86  session = LibUtilities::SessionReader::CreateInstance(argc, argv);
87 
88  // Create Field I/O object.
90  AllocateSharedPtr(session->GetComm());
91 
92  // Get some information from the session
93  sessionName = session->GetSessionName();
94  scheme = session->GetSolverInfo("TimeIntegrationMethod");
95  nSteps = session->GetParameter("NumSteps");
96  delta_t = session->GetParameter("TimeStep");
97  epsilon = session->GetParameter("epsilon");
98  lambda = 1.0/delta_t/epsilon;
99 
100  // Read the geometry and the expansion information
101  graph = SpatialDomains::MeshGraph::Read(session);
102 
103  // Set up the field
105  AllocateSharedPtr(session, graph, session->GetVariable(0));
106 
108  fields[0] = field->UpdatePhys();
109 
110  // Get coordinates of physical points
111  unsigned int nq = field->GetNpoints();
112  Array<OneD,NekDouble> x0(nq), x1(nq), x2(nq);
113  field->GetCoords(x0,x1,x2);
114 
115  // Evaluate initial condition
117  = session->GetFunction("InitialConditions", "u");
118  icond->Evaluate(x0,x1,x2,0.0,field->UpdatePhys());
119 }
120 
122 {
123  session->Finalise();
124 }
125 
127 {
129  CreateInstance(scheme);
130 
131  ode.DefineImplicitSolve(&Diffusion::DoImplicitSolve, this);
132 
133  // Initialise the scheme for actual time integration scheme
134  u = IntScheme->InitializeScheme(delta_t, fields, 0.0, ode);
135 
136  // Zero field coefficients for initial guess for linear solver.
137  Vmath::Zero(field->GetNcoeffs(), field->UpdateCoeffs(), 1);
138 
139  for (int n = 0; n < nSteps; ++n)
140  {
141  fields = IntScheme->TimeIntegrate(n, delta_t, u, ode);
142  }
143  Vmath::Vcopy(field->GetNpoints(), fields[0], 1, field->UpdatePhys(), 1);
144 
145  WriteSolution();
146  ExactSolution();
147 
148 }
149 
150 
152  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
153  Array<OneD, Array<OneD, NekDouble> >&outarray,
154  const NekDouble time,
155  const NekDouble lambda)
156 {
158  factors[StdRegions::eFactorLambda] = 1.0/lambda/epsilon;
159 
160  for (int i = 0; i < inarray.num_elements(); ++i)
161  {
162  // Multiply RHS by 1.0/timestep/lambda
163  Vmath::Smul(field->GetNpoints(), -factors[StdRegions::eFactorLambda],
164  inarray [i], 1,
165  outarray[i], 1);
166 
167  // Solve a system of equations with Helmholtz solver
168  field->HelmSolve(outarray[i],
169  field->UpdateCoeffs(),
170  NullFlagList, factors);
171 
172  // Transform to physical space and store in solution vector
173  field->BwdTrans (field->GetCoeffs(), outarray[i]);
174  }
175 }
176 
178 {
179  // Write solution to file
180  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
181  = field->GetFieldDefinitions();
182  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
183  for(int i = 0; i < FieldDef.size(); ++i)
184  {
185  FieldDef[i]->m_fields.push_back("u");
186  field->AppendFieldData(FieldDef[i], FieldData[i]);
187  }
188  fld->Write(session->GetSessionName() + ".fld", FieldDef, FieldData);
189 
190 }
191 
192 
194 {
195  unsigned int nq = field->GetNpoints();
196  Array<OneD,NekDouble> x0(nq), x1(nq), x2(nq);
197  field->GetCoords(x0,x1,x2);
198 
200  session->GetFunction("ExactSolution",0);
201 
202  if(ex_sol)
203  {
204  // evaluate exact solution
205  Array<OneD, NekDouble> exact(nq);
206  ex_sol->Evaluate(x0, x1, x2, (nSteps)*delta_t, exact);
207 
208  // Calculate errors
209  cout << "L inf error: "
210  << field->Linf(field->GetPhys(), exact) << endl;
211  cout << "L 2 error: "
212  << field->L2(field->GetPhys(), exact) << endl;
213  cout << "H 1 error: "
214  << field->H1(field->GetPhys(), exact) << endl;
215  }
216 
217 }
218 
219 int main(int argc, char *argv[])
220 {
221  try
222  {
223  Diffusion ops(argc, argv);
224  ops.TimeIntegrate();
225  }
226  catch (const std::runtime_error& e)
227  {
228  exit(-1);
229  }
230  catch (const std::string& eStr)
231  {
232  cout << "Error: " << eStr << endl;
233  exit(-1);
234  }
235 }
236 
LibUtilities::TimeIntegrationWrapperSharedPtr IntScheme
MultiRegions::ContField2DSharedPtr field
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.
Diffusion(int argc, char *argv[])
int main(int argc, char *argv[])
boost::shared_ptr< TimeIntegrationWrapper > TimeIntegrationWrapperSharedPtr
LibUtilities::TimeIntegrationSolutionSharedPtr u
LibUtilities::FieldIOSharedPtr fld
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
Array< OneD, Array< OneD, NekDouble > > fields
boost::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:236
double NekDouble
TimeIntegrationWrapperFactory & GetTimeIntegrationWrapperFactory()
boost::shared_ptr< Equation > EquationSharedPtr
LibUtilities::SessionReaderSharedPtr session
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble lambda)
SpatialDomains::MeshGraphSharedPtr graph
LibUtilities::TimeIntegrationSchemeOperators ode
unsigned int nSteps
boost::shared_ptr< TimeIntegrationSolution > TimeIntegrationSolutionSharedPtr
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
static FlagList NullFlagList
An empty flag list.