Nektar++
EigenValuesAdvection.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File EigenValuesAdvection.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:
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 
39 
40 namespace Nektar
41 {
42  string EigenValuesAdvection::className = GetEquationSystemFactory().RegisterCreatorFunction("EigenValuesAdvection", EigenValuesAdvection::create, "Eigenvalues of the weak advection operator.");
43 
46  : EquationSystem(pSession)
47  {
48  }
49 
51  {
52  EquationSystem::v_InitObject();
53 
54  // Define Velocity fields
56  std::vector<std::string> vel;
57  vel.push_back("Vx");
58  vel.push_back("Vy");
59  vel.push_back("Vz");
60  vel.resize(m_spacedim);
61 
62  EvaluateFunction(vel, m_velocity, "AdvectionVelocity");
63  }
64 
66  {
67 
68  }
69 
71  {
72 
73  }
74 
76  {
77  int nvariables = 1;
78  int i,dofs = GetNcoeffs();
79  //bool UseContCoeffs = false;
80 
81  Array<OneD, Array<OneD, NekDouble> > inarray(nvariables);
82  Array<OneD, Array<OneD, NekDouble> > tmp(nvariables);
83  Array<OneD, Array<OneD, NekDouble> > outarray(nvariables);
84  Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
85 
86  int npoints = GetNpoints();
87  int ncoeffs = GetNcoeffs();
88 
89  switch (m_projectionType)
90  {
92  {
93  dofs = ncoeffs;
94  break;
95  }
98  {
99  //dofs = GetContNcoeffs();
100  //UseContCoeffs = true;
101  break;
102  }
103  }
104 
105  cout << endl;
106  cout << "Num Phys Points = " << npoints << endl; // phisical points
107  cout << "Num Coeffs = " << ncoeffs << endl; //
108  cout << "Num Cont Coeffs = " << dofs << endl;
109 
110  inarray[0] = Array<OneD, NekDouble>(npoints,0.0);
111  outarray[0] = Array<OneD, NekDouble>(npoints,0.0);
112  tmp[0] = Array<OneD, NekDouble>(npoints,0.0);
113 
114  WeakAdv[0] = Array<OneD, NekDouble>(ncoeffs,0.0);
115  Array<OneD, NekDouble> MATRIX(npoints*npoints,0.0);
116 
117  for (int j = 0; j < npoints; j++)
118  {
119 
120  inarray[0][j] = 1.0;
121 
122  /// Feeding the weak Advection oprator with a vector (inarray)
123  /// Looping on inarray and changing the position of the only non-zero entry
124  /// we simulate the multiplication by the identity matrix.
125  /// The results stored in outarray is one of the columns of the weak advection oprators
126  /// which are then stored in MATRIX for the futher eigenvalues calculation.
127 
128  switch (m_projectionType)
129  {
131  {
132  WeakDGAdvection(inarray, WeakAdv,true,true,1);
133 
134  m_fields[0]->MultiplyByElmtInvMass(WeakAdv[0],WeakAdv[0]);
135 
136  m_fields[0]->BwdTrans(WeakAdv[0],outarray[0]);
137 
138  Vmath::Neg(npoints,outarray[0],1);
139  break;
140  }
143  {
144  // Calculate -V\cdot Grad(u);
145  for(i = 0; i < nvariables; ++i)
146  {
147  //Projection
148  m_fields[i]->FwdTrans(inarray[i],WeakAdv[i]);
149 
150  m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],tmp[i]);
151 
152  //Advection operator
153  AdvectionNonConservativeForm(m_velocity,tmp[i],outarray[i]);
154 
155  Vmath::Neg(npoints,outarray[i],1);
156 
157  //m_fields[i]->MultiplyByInvMassMatrix(WeakAdv[i],WeakAdv[i]);
158  //Projection
159  m_fields[i]->FwdTrans(outarray[i],WeakAdv[i]);
160 
161  m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],outarray[i]);
162  }
163  break;
164  }
165  }
166 
167  /// The result is stored in outarray (is the j-th columns of the weak advection operator).
168  /// We now store it in MATRIX(j)
169  Vmath::Vcopy(npoints,&(outarray[0][0]),1,&(MATRIX[j]),npoints);
170 
171  /// Set the j-th entry of inarray back to zero
172  inarray[0][j] = 0.0;
173  }
174 
175  ////////////////////////////////////////////////////////////////////////////////
176  /// Calulating the eigenvalues of the weak advection operator stored in (MATRIX)
177  /// using Lapack routines
178 
179  char jobvl = 'N';
180  char jobvr = 'N';
181  int info = 0, lwork = 3*npoints;
182  NekDouble dum;
183 
184  Array<OneD, NekDouble> EIG_R(npoints);
185  Array<OneD, NekDouble> EIG_I(npoints);
186 
187  Array<OneD, NekDouble> work(lwork);
188 
189  Lapack::Dgeev(jobvl,jobvr,npoints,MATRIX.get(),npoints,EIG_R.get(),EIG_I.get(),&dum,1,&dum,1,&work[0],lwork,info);
190 
191  ////////////////////////////////////////////////////////
192  //Print Matrix
193  FILE *mFile;
194 
195  mFile = fopen ("WeakAdvMatrix.txt","w");
196  for(int j = 0; j<npoints; j++)
197  {
198  for(int k = 0; k<npoints; k++)
199  {
200  fprintf(mFile,"%e ",MATRIX[j*npoints+k]);
201  }
202  fprintf(mFile,"\n");
203  }
204  fclose (mFile);
205 
206  ////////////////////////////////////////////////////////
207  //Output of the EigenValues
208  FILE *pFile;
209 
210  pFile = fopen ("Eigenvalues.txt","w");
211  for(int j = 0; j<npoints; j++)
212  {
213  fprintf(pFile,"%e %e\n",EIG_R[j],EIG_I[j]);
214  }
215  fclose (pFile);
216 
217  cout << "\nEigenvalues : " << endl;
218  for(int j = 0; j<npoints; j++)
219  {
220  cout << EIG_R[j] << "\t" << EIG_I[j] << endl;
221  }
222  cout << endl;
223  }
224 
227  {
228  ASSERTL1(flux.num_elements() == m_velocity.num_elements(),"Dimension of flux array and velocity array do not match");
229 
230  for(int j = 0; j < flux.num_elements(); ++j)
231  {
232  Vmath::Vmul(GetNpoints(),physfield[i],1,
233  m_velocity[j],1,flux[j],1);
234  }
235  }
236 
238  {
239  int i;
240 
241  int nTraceNumPoints = GetTraceNpoints();
242  int nvel = m_spacedim; //m_velocity.num_elements();
243 
244  Array<OneD, NekDouble > Fwd(nTraceNumPoints);
245  Array<OneD, NekDouble > Bwd(nTraceNumPoints);
246  Array<OneD, NekDouble > Vn (nTraceNumPoints,0.0);
247 
248  //Get Edge Velocity - Could be stored if time independent
249  for(i = 0; i < nvel; ++i)
250  {
251  m_fields[0]->ExtractTracePhys(m_velocity[i], Fwd);
252  Vmath::Vvtvp(nTraceNumPoints,m_traceNormals[i],1,Fwd,1,Vn,1,Vn,1);
253  }
254 
255  for(i = 0; i < numflux.num_elements(); ++i)
256  {
257  m_fields[i]->GetFwdBwdTracePhys(physfield[i],Fwd,Bwd);
258  m_fields[i]->GetTrace()->Upwind(Vn,Fwd,Bwd,numflux[i]);
259  // calculate m_fields[i]*Vn
260  Vmath::Vmul(nTraceNumPoints,numflux[i],1,Vn,1,numflux[i],1);
261  }
262  }
263 }
f SOLVER_UTILS_EXPORT void AdvectionNonConservativeForm(const Array< OneD, Array< OneD, NekDouble > > &V, const Array< OneD, const NekDouble > &u, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
Compute the non-conservative advection.
A base class for describing how to solve specific equations.
virtual void v_GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
static std::string className
Name of class.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
int m_spacedim
Spatial dimension (>= expansion dim).
EigenValuesAdvection(const LibUtilities::SessionReaderSharedPtr &pSession)
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
virtual void v_DoInitialise()
Virtual function for initialisation implementation.
virtual void v_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
SOLVER_UTILS_EXPORT void EvaluateFunction(Array< OneD, Array< OneD, NekDouble > > &pArray, std::string pFunctionName, const NekDouble pTime=0.0, const int domain=0)
Evaluates a function as specified in the session file.
EquationSystemFactory & GetEquationSystemFactory()
virtual void v_InitObject()
Initialisation object for EquationSystem.
Array< OneD, Array< OneD, NekDouble > > m_velocity
SOLVER_UTILS_EXPORT int GetNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SOLVER_UTILS_EXPORT int GetNcoeffs()
SOLVER_UTILS_EXPORT void WeakDGAdvection(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
Calculate the weak discontinuous Galerkin advection.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
virtual void v_DoSolve()
Virtual function for solve implementation.
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215