Nektar++
VelocityCorrectionScheme.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File VelocityCorrection.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: Velocity Correction Scheme for the Incompressible
33 // Navier Stokes equations
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 #include <SolverUtils/Core/Misc.h>
39 
40 #include <boost/algorithm/string.hpp>
41 
42 namespace Nektar
43 {
46  "VelocityCorrectionScheme",
48 
49  /**
50  * Constructor. Creates ...
51  *
52  * \param
53  * \param
54  */
57  : UnsteadySystem(pSession),
58  IncNavierStokes(pSession)
59  {
60 
61  }
62 
64  {
65  int n;
66 
68 
69  // Set m_pressure to point to last field of m_fields;
70  if (boost::iequals(m_session->GetVariable(m_fields.num_elements()-1), "p"))
71  {
72  m_nConvectiveFields = m_fields.num_elements()-1;
74  }
75  else
76  {
77  ASSERTL0(false,"Need to set up pressure field definition");
78  }
79 
80  // creation of the extrapolation object
82  {
83  std::string vExtrapolation = "Standard";
84 
85  if (m_session->DefinesSolverInfo("Extrapolation"))
86  {
87  vExtrapolation = m_session->GetSolverInfo("Extrapolation");
88  }
89 
91  vExtrapolation,
92  m_session,
93  m_fields,
94  m_pressure,
95  m_velocity,
96  m_advObject);
97  }
98 
99  // Integrate only the convective fields
100  for (n = 0; n < m_nConvectiveFields; ++n)
101  {
102  m_intVariables.push_back(n);
103  }
104 
105  // Load parameters for Spectral Vanishing Viscosity
106  m_session->MatchSolverInfo("SpectralVanishingViscosity","True",m_useSpecVanVisc,false);
107  m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
108  m_session->LoadParameter("SVVDiffCoeff",m_sVVDiffCoeff,0.1);
109  m_session->MatchSolverInfo("SPECTRALHPDEALIASING","True",m_specHP_dealiasing,false);
110 
111  // Needs to be set outside of next if so that it is turned off by default
112  m_session->MatchSolverInfo("SpectralVanishingViscosityHomo1D","True",m_useHomo1DSpecVanVisc,false);
113 
115  {
116  ASSERTL0(m_nConvectiveFields > 2,"Expect to have three velocity fields with homogenous expansion");
117 
118  if(m_useHomo1DSpecVanVisc == false)
119  {
120  m_session->MatchSolverInfo("SpectralVanishingViscosity","True",m_useHomo1DSpecVanVisc,false);
121  }
122 
124  {
126  planes = m_fields[0]->GetZIDs();
127 
128  int num_planes = planes.num_elements();
129  Array<OneD, NekDouble> SVV(num_planes,0.0);
130  NekDouble fac;
131  int kmodes = m_fields[0]->GetHomogeneousBasis()->GetNumModes();
132  int pstart;
133 
134  pstart = m_sVVCutoffRatio*kmodes;
135 
136  for(n = 0; n < num_planes; ++n)
137  {
138  if(planes[n] > pstart)
139  {
140  fac = (NekDouble)((planes[n] - kmodes)*(planes[n] - kmodes))/
141  ((NekDouble)((planes[n] - pstart)*(planes[n] - pstart)));
142  SVV[n] = m_sVVDiffCoeff*exp(-fac)/m_kinvis;
143  }
144  }
145 
146  for(int i = 0; i < m_velocity.num_elements(); ++i)
147  {
148  m_fields[m_velocity[i]]->SetHomo1DSpecVanVisc(SVV);
149  }
150  }
151 
152  }
153 
154  m_session->MatchSolverInfo("SmoothAdvection", "True", m_SmoothAdvection, false);
155 
156  if(m_subSteppingScheme) // Substepping
157  {
159  "Projection must be set to Mixed_CG_Discontinuous for "
160  "substepping");
161  }
162  else // Standard velocity correction scheme
163  {
164  // set explicit time-intregration class operators
166  }
167  m_extrapolation->SubSteppingTimeIntegration(m_intScheme->GetIntegrationMethod(), m_intScheme);
168  m_extrapolation->GenerateHOPBCMap();
169 
170  // set implicit time-intregration class operators
172  }
173 
174  /**
175  * Destructor
176  */
178  {
179  }
180 
181  /**
182  *
183  */
185  {
186  UnsteadySystem::v_GenerateSummary(s);
187 
189  {
191  s, "Substepping", LibUtilities::TimeIntegrationMethodMap[
192  m_subStepIntegrationScheme->GetIntegrationMethod()]);
193  }
194 
195  string dealias = m_homogen_dealiasing ? "Homogeneous1D" : "";
197  {
198  dealias += (dealias == "" ? "" : " + ") + string("spectral/hp");
199  }
200  if (dealias != "")
201  {
202  SolverUtils::AddSummaryItem(s, "Dealiasing", dealias);
203  }
204 
205  string smoothing = m_useSpecVanVisc ? "spectral/hp" : "";
207  {
208  smoothing += (smoothing == "" ? "" : " + ") + string("Homogeneous1D");
209  }
210  if (smoothing != "")
211  {
213  s, "Smoothing", "SVV (" + smoothing + " SVV (cut-off = "
214  + boost::lexical_cast<string>(m_sVVCutoffRatio)
215  + ", diff coeff = "
216  + boost::lexical_cast<string>(m_sVVDiffCoeff)+")");
217  }
218  }
219 
220  /**
221  *
222  */
224  {
225 
226  UnsteadySystem::v_DoInitialise();
227 
228  // Set up Field Meta Data for output files
229  m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
230  m_fieldMetaDataMap["TimeStep"] = boost::lexical_cast<std::string>(m_timestep);
231 
232  for(int i = 0; i < m_nConvectiveFields; ++i)
233  {
234  m_fields[i]->LocalToGlobal();
235  m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
236  m_fields[i]->GlobalToLocal();
237  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
238  m_fields[i]->UpdatePhys());
239  }
240  }
241 
242 
243  /**
244  *
245  */
247  {
248  int nfields = m_fields.num_elements() - 1;
249  for (int k=0 ; k < nfields; ++k)
250  {
251  //Backward Transformation in physical space for time evolution
252  m_fields[k]->BwdTrans_IterPerExp(m_fields[k]->GetCoeffs(),
253  m_fields[k]->UpdatePhys());
254  }
255  }
256 
257  /**
258  *
259  */
261  {
262 
263  int nfields = m_fields.num_elements() - 1;
264  for (int k=0 ; k < nfields; ++k)
265  {
266  //Forward Transformation in physical space for time evolution
267  m_fields[k]->FwdTrans_IterPerExp(m_fields[k]->GetPhys(),m_fields[k]->UpdateCoeffs());
268  }
269  }
270 
271  /**
272  *
273  */
275  {
276  int vVar = m_session->GetVariables().size();
277  Array<OneD, bool> vChecks(vVar, false);
278  vChecks[vVar-1] = true;
279  return vChecks;
280  }
281 
282  /**
283  *
284  */
286  {
287  return m_session->GetVariables().size() - 1;
288  }
289 
290  /**
291  * Explicit part of the method - Advection, Forcing + HOPBCs
292  */
294  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
295  Array<OneD, Array<OneD, NekDouble> > &outarray,
296  const NekDouble time)
297  {
298  EvaluateAdvectionTerms(inarray, outarray);
299 
300  // Smooth advection
302  {
303  for(int i = 0; i < m_nConvectiveFields; ++i)
304  {
305  m_pressure->SmoothField(outarray[i]);
306  }
307  }
308 
309  // Add forcing terms
310  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
311  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
312  {
313  (*x)->Apply(m_fields, inarray, outarray, time);
314  }
315 
316  // Calculate High-Order pressure boundary conditions
317  m_extrapolation->EvaluatePressureBCs(inarray,outarray,m_kinvis);
318  }
319 
320  /**
321  * Implicit part of the method - Poisson + nConv*Helmholtz
322  */
324  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
325  Array<OneD, Array<OneD, NekDouble> > &outarray,
326  const NekDouble time,
327  const NekDouble aii_Dt)
328  {
329  int i;
330  int phystot = m_fields[0]->GetTotPoints();
331 
333 
335  for(i = 0; i < m_nConvectiveFields; ++i)
336  {
337  F[i] = Array<OneD, NekDouble> (phystot);
338  }
339 
340  // Enforcing boundary conditions on all fields
341  SetBoundaryConditions(time);
342 
343  // Substep the pressure boundary condition
344  m_extrapolation->SubStepSetPressureBCs(inarray,aii_Dt,m_kinvis);
345 
346  // Set up forcing term and coefficients for pressure Poisson equation
347  SetUpPressureForcing(inarray, F, aii_Dt);
348  factors[StdRegions::eFactorLambda] = 0.0;
349 
350  // Solver Pressure Poisson Equation
351  m_pressure->HelmSolve(F[0], m_pressure->UpdateCoeffs(), NullFlagList,
352  factors);
353 
354  // Set up forcing term and coefficients for Helmholtz problems
355  SetUpViscousForcing(inarray, F, aii_Dt);
356  factors[StdRegions::eFactorLambda] = 1.0/aii_Dt/m_kinvis;
357  if(m_useSpecVanVisc)
358  {
361  }
362 
363  // Solve Helmholtz system and put in Physical space
364  for(i = 0; i < m_nConvectiveFields; ++i)
365  {
366  m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(),
367  NullFlagList, factors);
368  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),outarray[i]);
369  }
370  }
371 
372  /**
373  * Forcing term for Poisson solver solver
374  */
376  const Array<OneD, const Array<OneD, NekDouble> > &fields,
377  Array<OneD, Array<OneD, NekDouble> > &Forcing,
378  const NekDouble aii_Dt)
379  {
380  int i;
381  int physTot = m_fields[0]->GetTotPoints();
382  int nvel = m_velocity.num_elements();
383  Array<OneD, NekDouble> wk(physTot, 0.0);
384 
385  Vmath::Zero(physTot,Forcing[0],1);
386 
387  for(i = 0; i < nvel; ++i)
388  {
389  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[i],fields[i], wk);
390  Vmath::Vadd(physTot,wk,1,Forcing[0],1,Forcing[0],1);
391  }
392 
393  Vmath::Smul(physTot,1.0/aii_Dt,Forcing[0],1,Forcing[0],1);
394  }
395 
396  /**
397  * Forcing term for Helmholtz solver
398  */
400  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
401  Array<OneD, Array<OneD, NekDouble> > &Forcing,
402  const NekDouble aii_Dt)
403  {
404  NekDouble aii_dtinv = 1.0/aii_Dt;
405  int phystot = m_fields[0]->GetTotPoints();
406 
407  // Grad p
408  m_pressure->BwdTrans(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
409 
410  int nvel = m_velocity.num_elements();
411  if(nvel == 2)
412  {
413  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1]);
414  }
415  else
416  {
417  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1],
418  Forcing[2]);
419  }
420 
421  // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
422  // need to be updated for the convected fields.
423  for(int i = 0; i < nvel; ++i)
424  {
425  Blas::Daxpy(phystot,-aii_dtinv,inarray[i],1,Forcing[i],1);
426  Blas::Dscal(phystot,1.0/m_kinvis,&(Forcing[i])[0],1);
427  }
428  }
429 
430 } //end of namespace
EquationType m_equationType
equation type;
static SolverUtils::EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
virtual Array< OneD, bool > v_GetSystemSingularChecks()
virtual void v_DoInitialise(void)
Sets up initial conditions.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
void SetBoundaryConditions(NekDouble time)
time dependent boundary conditions updating
void EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
void SolveUnsteadyStokesSystem(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time, const NekDouble a_iixDt)
NekDouble m_kinvis
Kinematic viscosity.
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:48
VelocityCorrectionScheme(const LibUtilities::SessionReaderSharedPtr &pSession)
Constructor.
bool m_subSteppingScheme
bool to identify if using a substepping scheme
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
NekDouble m_timestep
Time step size.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
Array< OneD, int > m_velocity
int which identifies which components of m_fields contains the velocity (u,v,w);
bool m_useHomo1DSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
NekDouble m_sVVDiffCoeff
Diffusion coefficient of SVV modes.
ExtrapolateSharedPtr m_extrapolation
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
bool m_specHP_dealiasing
Flag to determine if dealisising is usde for the Spectral/hp element discretisation.
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
static std::string className
Name of class.
const char *const TimeIntegrationMethodMap[]
int m_nConvectiveFields
Number of fields to be convected;.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
void SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
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
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
std::vector< SolverUtils::ForcingSharedPtr > m_forcing
Forcing terms.
void SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
bool m_SmoothAdvection
bool to identify if advection term smoothing is requested
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:50
virtual void v_TransCoeffToPhys(void)
Virtual function for transformation to physical space.
LibUtilities::FieldMetaDataMap m_fieldMetaDataMap
Map to identify relevant solver info to dump in output fields.
double NekDouble
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
EquationSystemFactory & GetEquationSystemFactory()
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, Array< OneD, NekDouble > &wk=NullNekDouble1DArray)
NekDouble m_sVVCutoffRatio
cutt off ratio from which to start decayhing modes
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
bool m_useSpecVanVisc
bool to identify if spectral vanishing viscosity is active.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
virtual void v_TransPhysToCoeff(void)
Virtual function for transformation to coefficient space.
LibUtilities::TimeIntegrationWrapperSharedPtr m_intScheme
Wrapper to the time integration scheme.
This class is the base class for Navier Stokes problems.
LibUtilities::TimeIntegrationWrapperSharedPtr m_subStepIntegrationScheme
virtual void v_InitObject()
Init object for UnsteadySystem class.
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
virtual void v_InitObject()
Init object for UnsteadySystem class.
enum HomogeneousType m_HomogeneousType
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
static FlagList NullFlagList
An empty flag list.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215