Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 std;
45 using namespace Nektar;
46 
47 class Diffusion
48 {
49  public:
50  Diffusion(int argc, char* argv[]);
51  ~Diffusion();
52 
53  void TimeIntegrate();
54 
55  void DoImplicitSolve(
56  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
57  Array<OneD, Array<OneD, NekDouble> >&outarray,
58  const NekDouble time,
59  const NekDouble lambda);
60 
61  private:
64  string sessionName;
67 
72 
73  string scheme;
74  unsigned int nSteps;
78 
79  void WriteSolution();
80  void ExactSolution();
81 };
82 
83 
84 Diffusion::Diffusion(int argc, char* argv[])
85 {
86  // Create session reader.
87  session = LibUtilities::SessionReader::CreateInstance(argc, argv);
88 
89  // Create Field I/O object.
91  AllocateSharedPtr(session->GetComm());
92 
93  // Get some information from the session
94  sessionName = session->GetSessionName();
95  scheme = session->GetSolverInfo("TimeIntegrationMethod");
96  nSteps = session->GetParameter("NumSteps");
97  delta_t = session->GetParameter("TimeStep");
98  epsilon = session->GetParameter("epsilon");
99  lambda = 1.0/delta_t/epsilon;
100 
101  // Read the geometry and the expansion information
102  graph = SpatialDomains::MeshGraph::Read(session);
103 
104  // Set up the field
106  AllocateSharedPtr(session, graph, session->GetVariable(0));
107 
109  fields[0] = field->UpdatePhys();
110 
111  // Get coordinates of physical points
112  unsigned int nq = field->GetNpoints();
113  Array<OneD,NekDouble> x0(nq), x1(nq), x2(nq);
114  field->GetCoords(x0,x1,x2);
115 
116  // Evaluate initial condition
118  = session->GetFunction("InitialConditions", "u");
119  icond->Evaluate(x0,x1,x2,0.0,field->UpdatePhys());
120 }
121 
123 {
124  session->Finalise();
125 }
126 
128 {
130  CreateInstance(scheme);
131 
132  ode.DefineImplicitSolve(&Diffusion::DoImplicitSolve, this);
133 
134  // Initialise the scheme for actual time integration scheme
135  u = IntScheme->InitializeScheme(delta_t, fields, 0.0, ode);
136 
137  // Zero field coefficients for initial guess for linear solver.
138  Vmath::Zero(field->GetNcoeffs(), field->UpdateCoeffs(), 1);
139 
140  for (int n = 0; n < nSteps; ++n)
141  {
142  fields = IntScheme->TimeIntegrate(n, delta_t, u, ode);
143  }
144  Vmath::Vcopy(field->GetNpoints(), fields[0], 1, field->UpdatePhys(), 1);
145 
146  WriteSolution();
147  ExactSolution();
148 
149 }
150 
151 
153  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
154  Array<OneD, Array<OneD, NekDouble> >&outarray,
155  const NekDouble time,
156  const NekDouble lambda)
157 {
159  factors[StdRegions::eFactorLambda] = 1.0/lambda/epsilon;
160 
161  for (int i = 0; i < inarray.num_elements(); ++i)
162  {
163  // Multiply RHS by 1.0/timestep/lambda
164  Vmath::Smul(field->GetNpoints(), -factors[StdRegions::eFactorLambda],
165  inarray [i], 1,
166  outarray[i], 1);
167 
168  // Solve a system of equations with Helmholtz solver
169  field->HelmSolve(outarray[i],
170  field->UpdateCoeffs(),
171  NullFlagList, factors);
172 
173  // Transform to physical space and store in solution vector
174  field->BwdTrans (field->GetCoeffs(), outarray[i]);
175  }
176 }
177 
179 {
180  // Write solution to file
181  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
182  = field->GetFieldDefinitions();
183  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
184  for(int i = 0; i < FieldDef.size(); ++i)
185  {
186  FieldDef[i]->m_fields.push_back("u");
187  field->AppendFieldData(FieldDef[i], FieldData[i]);
188  }
189  fld->Write(session->GetSessionName() + ".fld", FieldDef, FieldData);
190 
191 }
192 
193 
195 {
196  unsigned int nq = field->GetNpoints();
197  Array<OneD,NekDouble> x0(nq), x1(nq), x2(nq);
198  field->GetCoords(x0,x1,x2);
199 
201  session->GetFunction("ExactSolution",0);
202 
203  if(ex_sol)
204  {
205  // evaluate exact solution
206  Array<OneD, NekDouble> exact(nq);
207  ex_sol->Evaluate(x0, x1, x2, (nSteps)*delta_t, exact);
208 
209  // Calculate errors
210  cout << "L inf error: "
211  << field->Linf(field->GetPhys(), exact) << endl;
212  cout << "L 2 error: "
213  << field->L2(field->GetPhys(), exact) << endl;
214  cout << "H 1 error: "
215  << field->H1(field->GetPhys(), exact) << endl;
216  }
217 
218 }
219 
220 int main(int argc, char *argv[])
221 {
222  try
223  {
224  Diffusion ops(argc, argv);
225  ops.TimeIntegrate();
226  }
227  catch (const std::runtime_error& e)
228  {
229  exit(-1);
230  }
231  catch (const std::string& eStr)
232  {
233  cout << "Error: " << eStr << endl;
234  exit(-1);
235  }
236 }
237 
LibUtilities::TimeIntegrationWrapperSharedPtr IntScheme
MultiRegions::ContField2DSharedPtr field
Diffusion(int argc, char *argv[])
int main(int argc, char *argv[])
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
boost::shared_ptr< TimeIntegrationWrapper > TimeIntegrationWrapperSharedPtr
LibUtilities::TimeIntegrationSolutionSharedPtr u
STL namespace.
LibUtilities::FieldIOSharedPtr fld
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:293
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:251
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
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:225
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:442
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
static FlagList NullFlagList
An empty flag list.