Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 using namespace std;
43 
44 namespace Nektar
45 {
46  string VelocityCorrectionScheme::className =
48  "VelocityCorrectionScheme",
49  VelocityCorrectionScheme::create);
50 
51  /**
52  * Constructor. Creates ...
53  *
54  * \param
55  * \param
56  */
57  VelocityCorrectionScheme::VelocityCorrectionScheme(
59  : UnsteadySystem(pSession),
60  IncNavierStokes(pSession),
61  m_varCoeffLap(StdRegions::NullVarCoeffMap)
62  {
63 
64  }
65 
67  {
68  int n;
69 
71  m_explicitDiffusion = false;
72 
73  // Set m_pressure to point to last field of m_fields;
74  if (boost::iequals(m_session->GetVariable(m_fields.num_elements()-1), "p"))
75  {
76  m_nConvectiveFields = m_fields.num_elements()-1;
78  }
79  else
80  {
81  ASSERTL0(false,"Need to set up pressure field definition");
82  }
83 
84  // creation of the extrapolation object
86  {
87  std::string vExtrapolation = "Standard";
88 
89  if (m_session->DefinesSolverInfo("Extrapolation"))
90  {
91  vExtrapolation = m_session->GetSolverInfo("Extrapolation");
92  }
93 
95  vExtrapolation,
96  m_session,
97  m_fields,
98  m_pressure,
99  m_velocity,
100  m_advObject);
101  }
102 
103  // Integrate only the convective fields
104  for (n = 0; n < m_nConvectiveFields; ++n)
105  {
106  m_intVariables.push_back(n);
107  }
108 
111 
112  // Load parameters for Spectral Vanishing Viscosity
113  m_session->MatchSolverInfo("SpectralVanishingViscosity","True",
114  m_useSpecVanVisc,false);
115  m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
116  m_session->LoadParameter("SVVDiffCoeff", m_sVVDiffCoeff, 0.1);
117  // Needs to be set outside of next if so that it is turned off by default
118  m_session->MatchSolverInfo("SpectralVanishingViscosityHomo1D","True",
119  m_useHomo1DSpecVanVisc,false);
120 
121 
122  m_session->MatchSolverInfo("SPECTRALHPDEALIASING","True",
123  m_specHP_dealiasing,false);
124 
125 
127  {
128  ASSERTL0(m_nConvectiveFields > 2,"Expect to have three velocity fields with homogenous expansion");
129 
130  if(m_useHomo1DSpecVanVisc == false)
131  {
132  m_session->MatchSolverInfo("SpectralVanishingViscosity","True",m_useHomo1DSpecVanVisc,false);
133  }
134 
136  {
138  planes = m_fields[0]->GetZIDs();
139 
140  int num_planes = planes.num_elements();
141  Array<OneD, NekDouble> SVV(num_planes,0.0);
142  NekDouble fac;
143  int kmodes = m_fields[0]->GetHomogeneousBasis()->GetNumModes();
144  int pstart;
145 
146  pstart = m_sVVCutoffRatio*kmodes;
147 
148  for(n = 0; n < num_planes; ++n)
149  {
150  if(planes[n] > pstart)
151  {
152  fac = (NekDouble)((planes[n] - kmodes)*(planes[n] - kmodes))/
153  ((NekDouble)((planes[n] - pstart)*(planes[n] - pstart)));
154  SVV[n] = m_sVVDiffCoeff*exp(-fac)/m_kinvis;
155  }
156  }
157 
158  for(int i = 0; i < m_velocity.num_elements(); ++i)
159  {
160  m_fields[m_velocity[i]]->SetHomo1DSpecVanVisc(SVV);
161  }
162  }
163 
164  }
165 
166  m_session->MatchSolverInfo("SmoothAdvection", "True", m_SmoothAdvection, false);
167 
168  // set explicit time-intregration class operators
170 
171  m_extrapolation->SubSteppingTimeIntegration(m_intScheme->GetIntegrationMethod(), m_intScheme);
172  m_extrapolation->GenerateHOPBCMap();
173 
174  // set implicit time-intregration class operators
176  }
177 
178  /**
179  * Destructor
180  */
182  {
183  }
184 
185  /**
186  *
187  */
189  {
191 
192  if (m_extrapolation->GetSubStepIntegrationMethod() !=
194  {
195  SolverUtils::AddSummaryItem(s, "Substepping",
197  m_extrapolation->GetSubStepIntegrationMethod()]);
198  }
199 
200  string dealias = m_homogen_dealiasing ? "Homogeneous1D" : "";
202  {
203  dealias += (dealias == "" ? "" : " + ") + string("spectral/hp");
204  }
205  if (dealias != "")
206  {
207  SolverUtils::AddSummaryItem(s, "Dealiasing", dealias);
208  }
209 
210  string smoothing = m_useSpecVanVisc ? "spectral/hp" : "";
212  {
213  smoothing += (smoothing == "" ? "" : " + ") + string("Homogeneous1D");
214  }
215  if (smoothing != "")
216  {
218  s, "Smoothing", "SVV (" + smoothing + " SVV (cut-off = "
219  + boost::lexical_cast<string>(m_sVVCutoffRatio)
220  + ", diff coeff = "
221  + boost::lexical_cast<string>(m_sVVDiffCoeff)+")");
222  }
223  }
224 
225  /**
226  *
227  */
229  {
230 
232 
233  // Set up Field Meta Data for output files
234  m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
235  m_fieldMetaDataMap["TimeStep"] = boost::lexical_cast<std::string>(m_timestep);
236 
237  // set boundary conditions here so that any normal component
238  // correction are imposed before they are imposed on intiial
239  // field below
241 
242  for(int i = 0; i < m_nConvectiveFields; ++i)
243  {
244  m_fields[i]->LocalToGlobal();
245  m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
246  m_fields[i]->GlobalToLocal();
247  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
248  m_fields[i]->UpdatePhys());
249  }
250  }
251 
252 
253  /**
254  *
255  */
257  {
258  int nfields = m_fields.num_elements() - 1;
259  for (int k=0 ; k < nfields; ++k)
260  {
261  //Backward Transformation in physical space for time evolution
262  m_fields[k]->BwdTrans_IterPerExp(m_fields[k]->GetCoeffs(),
263  m_fields[k]->UpdatePhys());
264  }
265  }
266 
267  /**
268  *
269  */
271  {
272 
273  int nfields = m_fields.num_elements() - 1;
274  for (int k=0 ; k < nfields; ++k)
275  {
276  //Forward Transformation in physical space for time evolution
277  m_fields[k]->FwdTrans_IterPerExp(m_fields[k]->GetPhys(),m_fields[k]->UpdateCoeffs());
278  }
279  }
280 
281  /**
282  *
283  */
285  {
286  int vVar = m_session->GetVariables().size();
287  Array<OneD, bool> vChecks(vVar, false);
288  vChecks[vVar-1] = true;
289  return vChecks;
290  }
291 
292  /**
293  *
294  */
296  {
297  return m_session->GetVariables().size() - 1;
298  }
299 
300  /**
301  * Explicit part of the method - Advection, Forcing + HOPBCs
302  */
304  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
305  Array<OneD, Array<OneD, NekDouble> > &outarray,
306  const NekDouble time)
307  {
308  EvaluateAdvectionTerms(inarray, outarray);
309 
310  // Smooth advection
312  {
313  for(int i = 0; i < m_nConvectiveFields; ++i)
314  {
315  m_pressure->SmoothField(outarray[i]);
316  }
317  }
318 
319  // Add forcing terms
320  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
321  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
322  {
323  (*x)->Apply(m_fields, inarray, outarray, time);
324  }
325 
326  // Calculate High-Order pressure boundary conditions
327  m_extrapolation->EvaluatePressureBCs(inarray,outarray,m_kinvis);
328  }
329 
330  /**
331  * Implicit part of the method - Poisson + nConv*Helmholtz
332  */
334  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
335  Array<OneD, Array<OneD, NekDouble> > &outarray,
336  const NekDouble time,
337  const NekDouble aii_Dt)
338  {
339  int i;
340  int phystot = m_fields[0]->GetTotPoints();
341 
343  for(i = 0; i < m_nConvectiveFields; ++i)
344  {
345  F[i] = Array<OneD, NekDouble> (phystot);
346  }
347 
348  // Enforcing boundary conditions on all fields
349  SetBoundaryConditions(time);
350 
351  // Substep the pressure boundary condition if using substepping
352  m_extrapolation->SubStepSetPressureBCs(inarray,aii_Dt,m_kinvis);
353 
354  // Set up forcing term for pressure Poisson equation
355  SetUpPressureForcing(inarray, F, aii_Dt);
356 
357  // Solve Pressure System
358  SolvePressure (F[0]);
359 
360  // Set up forcing term for Helmholtz problems
361  SetUpViscousForcing(inarray, F, aii_Dt);
362 
363  // Solve velocity system
364  SolveViscous( F, outarray, aii_Dt);
365  }
366 
367  /**
368  * Forcing term for Poisson solver solver
369  */
371  const Array<OneD, const Array<OneD, NekDouble> > &fields,
373  const NekDouble aii_Dt)
374  {
375  int i;
376  int physTot = m_fields[0]->GetTotPoints();
377  int nvel = m_velocity.num_elements();
378  Array<OneD, NekDouble> wk(physTot, 0.0);
379 
380  Vmath::Zero(physTot,Forcing[0],1);
381 
382  for(i = 0; i < nvel; ++i)
383  {
384  m_fields[i]->PhysDeriv(MultiRegions::DirCartesianMap[i],fields[i], wk);
385  Vmath::Vadd(physTot,wk,1,Forcing[0],1,Forcing[0],1);
386  }
387 
388  Vmath::Smul(physTot,1.0/aii_Dt,Forcing[0],1,Forcing[0],1);
389  }
390 
391  /**
392  * Forcing term for Helmholtz solver
393  */
395  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
397  const NekDouble aii_Dt)
398  {
399  NekDouble aii_dtinv = 1.0/aii_Dt;
400  int phystot = m_fields[0]->GetTotPoints();
401 
402  // Grad p
403  m_pressure->BwdTrans(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
404 
405  int nvel = m_velocity.num_elements();
406  if(nvel == 2)
407  {
408  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1]);
409  }
410  else
411  {
412  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0],
413  Forcing[1], Forcing[2]);
414  }
415 
416  // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
417  // need to be updated for the convected fields.
418  for(int i = 0; i < nvel; ++i)
419  {
420  Blas::Daxpy(phystot,-aii_dtinv,inarray[i],1,Forcing[i],1);
421  Blas::Dscal(phystot,1.0/m_kinvis,&(Forcing[i])[0],1);
422  }
423  }
424 
425 
426  /**
427  * Solve pressure system
428  */
431  {
433  // Setup coefficient for equation
434  factors[StdRegions::eFactorLambda] = 0.0;
435 
436  // Solver Pressure Poisson Equation
437  m_pressure->HelmSolve(Forcing, m_pressure->UpdateCoeffs(), NullFlagList,
438  factors);
439  }
440 
441  /**
442  * Solve velocity system
443  */
445  const Array<OneD, const Array<OneD, NekDouble> > &Forcing,
446  Array<OneD, Array<OneD, NekDouble> > &outarray,
447  const NekDouble aii_Dt)
448  {
450  // Setup coefficients for equation
451  factors[StdRegions::eFactorLambda] = 1.0/aii_Dt/m_kinvis;
452  if(m_useSpecVanVisc)
453  {
456  }
457 
458  // Solve Helmholtz system and put in Physical space
459  for(int i = 0; i < m_nConvectiveFields; ++i)
460  {
461  m_fields[i]->HelmSolve(Forcing[i], m_fields[i]->UpdateCoeffs(),
462  NullFlagList, factors);
463  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),outarray[i]);
464  }
465  }
466 
467 } //end of namespace
EquationType m_equationType
equation type;
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_time
Current time of simulation.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
NekDouble m_kinvis
Kinematic viscosity.
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
void SolvePressure(const Array< OneD, NekDouble > &Forcing)
ExtrapolateFactory & GetExtrapolateFactory()
Definition: Extrapolate.cpp:50
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
STL namespace.
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:251
const char *const TimeIntegrationMethodMap[]
int m_nConvectiveFields
Number of fields to be convected;.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
void SolveViscous(const Array< OneD, const Array< OneD, NekDouble > > &Forcing, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble aii_Dt)
void SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
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.
Base class for unsteady solvers.
virtual void v_SetUpViscousForcing(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
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
virtual void v_SolvePressure(const Array< OneD, NekDouble > &Forcing)
static const NekDouble kNekUnsetDouble
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
EquationSystemFactory & GetEquationSystemFactory()
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:86
virtual void v_SolveViscous(const Array< OneD, const Array< OneD, NekDouble > > &Forcing, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble aii_Dt)
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.
Array< OneD, NekDouble > m_saved_aii_Dt
Save aiiDt value to use as a trip to reset global matrix setup.
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.
virtual void v_InitObject()
Init object for UnsteadySystem class.
MultiRegions::ExpListSharedPtr m_pressure
Pointer to field holding pressure field.
Defines a forcing term to be explicitly applied.
Definition: Forcing.h:70
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
virtual void v_SetUpPressureForcing(const Array< OneD, const Array< OneD, NekDouble > > &fields, Array< OneD, Array< OneD, NekDouble > > &Forcing, const NekDouble aii_Dt)
virtual void v_InitObject()
Init object for UnsteadySystem class.
virtual void v_EvaluateAdvection_SetPressureBCs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
virtual SOLVER_UTILS_EXPORT void v_DoInitialise()
Sets up initial conditions.
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.
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:227
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215