Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 using namespace std;
41 
42 namespace Nektar
43 {
44  string EigenValuesAdvection::className = GetEquationSystemFactory().RegisterCreatorFunction("EigenValuesAdvection", EigenValuesAdvection::create, "Eigenvalues of the weak advection operator.");
45 
46  EigenValuesAdvection::EigenValuesAdvection(
48  : EquationSystem(pSession)
49  {
50  }
51 
53  {
55 
56  // Define Velocity fields
58  std::vector<std::string> vel;
59  vel.push_back("Vx");
60  vel.push_back("Vy");
61  vel.push_back("Vz");
62  vel.resize(m_spacedim);
63 
64  EvaluateFunction(vel, m_velocity, "AdvectionVelocity");
65  }
66 
68  {
69 
70  }
71 
73  {
74 
75  }
76 
78  {
79  int nvariables = 1;
80  int i,dofs = GetNcoeffs();
81  //bool UseContCoeffs = false;
82 
83  Array<OneD, Array<OneD, NekDouble> > inarray(nvariables);
84  Array<OneD, Array<OneD, NekDouble> > tmp(nvariables);
85  Array<OneD, Array<OneD, NekDouble> > outarray(nvariables);
86  Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
87 
88  int npoints = GetNpoints();
89  int ncoeffs = GetNcoeffs();
90 
91  switch (m_projectionType)
92  {
94  {
95  dofs = ncoeffs;
96  break;
97  }
100  {
101  //dofs = GetContNcoeffs();
102  //UseContCoeffs = true;
103  break;
104  }
105  }
106 
107  cout << endl;
108  cout << "Num Phys Points = " << npoints << endl; // phisical points
109  cout << "Num Coeffs = " << ncoeffs << endl; //
110  cout << "Num Cont Coeffs = " << dofs << endl;
111 
112  inarray[0] = Array<OneD, NekDouble>(npoints,0.0);
113  outarray[0] = Array<OneD, NekDouble>(npoints,0.0);
114  tmp[0] = Array<OneD, NekDouble>(npoints,0.0);
115 
116  WeakAdv[0] = Array<OneD, NekDouble>(ncoeffs,0.0);
117  Array<OneD, NekDouble> MATRIX(npoints*npoints,0.0);
118 
119  for (int j = 0; j < npoints; j++)
120  {
121 
122  inarray[0][j] = 1.0;
123 
124  /// Feeding the weak Advection oprator with a vector (inarray)
125  /// Looping on inarray and changing the position of the only non-zero entry
126  /// we simulate the multiplication by the identity matrix.
127  /// The results stored in outarray is one of the columns of the weak advection oprators
128  /// which are then stored in MATRIX for the futher eigenvalues calculation.
129 
130  switch (m_projectionType)
131  {
133  {
134  WeakDGAdvection(inarray, WeakAdv,true,true,1);
135 
136  m_fields[0]->MultiplyByElmtInvMass(WeakAdv[0],WeakAdv[0]);
137 
138  m_fields[0]->BwdTrans(WeakAdv[0],outarray[0]);
139 
140  Vmath::Neg(npoints,outarray[0],1);
141  break;
142  }
145  {
146  // Calculate -V\cdot Grad(u);
147  for(i = 0; i < nvariables; ++i)
148  {
149  //Projection
150  m_fields[i]->FwdTrans(inarray[i],WeakAdv[i]);
151 
152  m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],tmp[i]);
153 
154  //Advection operator
155  AdvectionNonConservativeForm(m_velocity,tmp[i],outarray[i]);
156 
157  Vmath::Neg(npoints,outarray[i],1);
158 
159  //m_fields[i]->MultiplyByInvMassMatrix(WeakAdv[i],WeakAdv[i]);
160  //Projection
161  m_fields[i]->FwdTrans(outarray[i],WeakAdv[i]);
162 
163  m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],outarray[i]);
164  }
165  break;
166  }
167  }
168 
169  /// The result is stored in outarray (is the j-th columns of the weak advection operator).
170  /// We now store it in MATRIX(j)
171  Vmath::Vcopy(npoints,&(outarray[0][0]),1,&(MATRIX[j]),npoints);
172 
173  /// Set the j-th entry of inarray back to zero
174  inarray[0][j] = 0.0;
175  }
176 
177  ////////////////////////////////////////////////////////////////////////////////
178  /// Calulating the eigenvalues of the weak advection operator stored in (MATRIX)
179  /// using Lapack routines
180 
181  char jobvl = 'N';
182  char jobvr = 'N';
183  int info = 0, lwork = 3*npoints;
184  NekDouble dum;
185 
186  Array<OneD, NekDouble> EIG_R(npoints);
187  Array<OneD, NekDouble> EIG_I(npoints);
188 
189  Array<OneD, NekDouble> work(lwork);
190 
191  Lapack::Dgeev(jobvl,jobvr,npoints,MATRIX.get(),npoints,EIG_R.get(),EIG_I.get(),&dum,1,&dum,1,&work[0],lwork,info);
192 
193  ////////////////////////////////////////////////////////
194  //Print Matrix
195  FILE *mFile;
196 
197  mFile = fopen ("WeakAdvMatrix.txt","w");
198  for(int j = 0; j<npoints; j++)
199  {
200  for(int k = 0; k<npoints; k++)
201  {
202  fprintf(mFile,"%e ",MATRIX[j*npoints+k]);
203  }
204  fprintf(mFile,"\n");
205  }
206  fclose (mFile);
207 
208  ////////////////////////////////////////////////////////
209  //Output of the EigenValues
210  FILE *pFile;
211 
212  pFile = fopen ("Eigenvalues.txt","w");
213  for(int j = 0; j<npoints; j++)
214  {
215  fprintf(pFile,"%e %e\n",EIG_R[j],EIG_I[j]);
216  }
217  fclose (pFile);
218 
219  cout << "\nEigenvalues : " << endl;
220  for(int j = 0; j<npoints; j++)
221  {
222  cout << EIG_R[j] << "\t" << EIG_I[j] << endl;
223  }
224  cout << endl;
225  }
226 
229  {
230  ASSERTL1(flux.num_elements() == m_velocity.num_elements(),"Dimension of flux array and velocity array do not match");
231 
232  for(int j = 0; j < flux.num_elements(); ++j)
233  {
234  Vmath::Vmul(GetNpoints(),physfield[i],1,
235  m_velocity[j],1,flux[j],1);
236  }
237  }
238 
240  {
241  int i;
242 
243  int nTraceNumPoints = GetTraceNpoints();
244  int nvel = m_spacedim; //m_velocity.num_elements();
245 
246  Array<OneD, NekDouble > Fwd(nTraceNumPoints);
247  Array<OneD, NekDouble > Bwd(nTraceNumPoints);
248  Array<OneD, NekDouble > Vn (nTraceNumPoints,0.0);
249 
250  //Get Edge Velocity - Could be stored if time independent
251  for(i = 0; i < nvel; ++i)
252  {
253  m_fields[0]->ExtractTracePhys(m_velocity[i], Fwd);
254  Vmath::Vvtvp(nTraceNumPoints,m_traceNormals[i],1,Fwd,1,Vn,1,Vn,1);
255  }
256 
257  for(i = 0; i < numflux.num_elements(); ++i)
258  {
259  m_fields[i]->GetFwdBwdTracePhys(physfield[i],Fwd,Bwd);
260  m_fields[i]->GetTrace()->Upwind(Vn,Fwd,Bwd,numflux[i]);
261  // calculate m_fields[i]*Vn
262  Vmath::Vmul(nTraceNumPoints,numflux[i],1,Vn,1,numflux[i],1);
263  }
264  }
265 }
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:442
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
STL namespace.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Initialisation object for EquationSystem.
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).
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
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:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
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:183
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215