Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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  using namespace MultiRegions;
47 
48  string VelocityCorrectionScheme::className =
50  "VelocityCorrectionScheme",
51  VelocityCorrectionScheme::create);
52 
53  /**
54  * Constructor. Creates ...
55  *
56  * \param
57  * \param
58  */
59  VelocityCorrectionScheme::VelocityCorrectionScheme(
61  : UnsteadySystem(pSession),
62  IncNavierStokes(pSession),
63  m_varCoeffLap(StdRegions::NullVarCoeffMap)
64  {
65 
66  }
67 
69  {
70  int n;
71 
73  m_explicitDiffusion = false;
74 
75  // Set m_pressure to point to last field of m_fields;
76  if (boost::iequals(m_session->GetVariable(m_fields.num_elements()-1), "p"))
77  {
78  m_nConvectiveFields = m_fields.num_elements()-1;
80  }
81  else
82  {
83  ASSERTL0(false,"Need to set up pressure field definition");
84  }
85 
86  // Determine diffusion coefficients for each field
88  for (n = 0; n < m_nConvectiveFields; ++n)
89  {
90  std::string varName = m_session->GetVariable(n);
91  if ( m_session->DefinesFunction("DiffusionCoefficient", varName))
92  {
94  = m_session->GetFunction("DiffusionCoefficient", varName);
95  m_diffCoeff[n] = ffunc->Evaluate();
96  }
97  }
98 
99  // creation of the extrapolation object
101  {
102  std::string vExtrapolation = v_GetExtrapolateStr();
103 
104  if (m_session->DefinesSolverInfo("Extrapolation"))
105  {
106  vExtrapolation = v_GetSubSteppingExtrapolateStr(
107  m_session->GetSolverInfo("Extrapolation"));
108  }
109 
111  vExtrapolation,
112  m_session,
113  m_fields,
114  m_pressure,
115  m_velocity,
116  m_advObject);
117  }
118 
119  // Integrate only the convective fields
120  for (n = 0; n < m_nConvectiveFields; ++n)
121  {
122  m_intVariables.push_back(n);
123  }
124 
125  // Load parameters for Spectral Vanishing Viscosity
126  m_session->MatchSolverInfo("SpectralVanishingViscosity","True",
127  m_useSpecVanVisc, false);
129  if(m_useSpecVanVisc == false)
130  {
131  m_session->MatchSolverInfo("SpectralVanishingViscositySpectralHP",
132  "True", m_useSpecVanVisc, false);
133  m_session->MatchSolverInfo("SpectralVanishingViscosityHomo1D",
134  "True", m_useHomo1DSpecVanVisc, false);
135  }
136  m_session->LoadParameter("SVVCutoffRatio",m_sVVCutoffRatio,0.75);
137  m_session->LoadParameter("SVVDiffCoeff", m_sVVDiffCoeff, 0.1);
138 
139  m_session->MatchSolverInfo("SPECTRALHPDEALIASING","True",
140  m_specHP_dealiasing,false);
141 
143  {
144  ASSERTL0(m_nConvectiveFields > 2,"Expect to have three velocity fields with homogenous expansion");
145 
147  {
149  planes = m_fields[0]->GetZIDs();
150 
151  int num_planes = planes.num_elements();
152  Array<OneD, NekDouble> SVV(num_planes,0.0);
153  NekDouble fac;
154  int kmodes = m_fields[0]->GetHomogeneousBasis()->GetNumModes();
155  int pstart;
156 
157  pstart = m_sVVCutoffRatio*kmodes;
158 
159  for(n = 0; n < num_planes; ++n)
160  {
161  if(planes[n] > pstart)
162  {
163  fac = (NekDouble)((planes[n] - kmodes)*(planes[n] - kmodes))/
164  ((NekDouble)((planes[n] - pstart)*(planes[n] - pstart)));
165  SVV[n] = m_sVVDiffCoeff*exp(-fac)/m_kinvis;
166  }
167  }
168 
169  for(int i = 0; i < m_velocity.num_elements(); ++i)
170  {
171  m_fields[m_velocity[i]]->SetHomo1DSpecVanVisc(SVV);
172  }
173  }
174 
175  }
176 
177  m_session->MatchSolverInfo("SmoothAdvection", "True", m_SmoothAdvection, false);
178 
179  // set explicit time-intregration class operators
181 
182  m_extrapolation->SubSteppingTimeIntegration(m_intScheme->GetIntegrationMethod(), m_intScheme);
183  m_extrapolation->GenerateHOPBCMap(m_session);
184 
185  // set implicit time-intregration class operators
187  }
188 
189  /**
190  * Destructor
191  */
193  {
194  }
195 
196  /**
197  *
198  */
200  {
202  SolverUtils::AddSummaryItem(s, "Splitting Scheme", "Velocity correction (strong press. form)");
203 
204  if (m_extrapolation->GetSubStepIntegrationMethod() !=
206  {
207  SolverUtils::AddSummaryItem(s, "Substepping",
209  m_extrapolation->GetSubStepIntegrationMethod()]);
210  }
211 
212  string dealias = m_homogen_dealiasing ? "Homogeneous1D" : "";
214  {
215  dealias += (dealias == "" ? "" : " + ") + string("spectral/hp");
216  }
217  if (dealias != "")
218  {
219  SolverUtils::AddSummaryItem(s, "Dealiasing", dealias);
220  }
221 
222  string smoothing = m_useSpecVanVisc ? "spectral/hp" : "";
224  {
225  smoothing += (smoothing == "" ? "" : " + ") + string("Homogeneous1D");
226  }
227  if (smoothing != "")
228  {
230  s, "Smoothing", "SVV (" + smoothing + " SVV (cut-off = "
231  + boost::lexical_cast<string>(m_sVVCutoffRatio)
232  + ", diff coeff = "
233  + boost::lexical_cast<string>(m_sVVDiffCoeff)+")");
234  }
235  }
236 
237  /**
238  *
239  */
241  {
242 
244 
245  // Set up Field Meta Data for output files
246  m_fieldMetaDataMap["Kinvis"] = boost::lexical_cast<std::string>(m_kinvis);
247  m_fieldMetaDataMap["TimeStep"] = boost::lexical_cast<std::string>(m_timestep);
248 
249  // set boundary conditions here so that any normal component
250  // correction are imposed before they are imposed on intiial
251  // field below
253 
255  for(int i = 0; i < m_nConvectiveFields; ++i)
256  {
257  m_fields[i]->LocalToGlobal();
258  m_fields[i]->ImposeDirichletConditions(m_fields[i]->UpdateCoeffs());
259  m_fields[i]->GlobalToLocal();
260  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
261  m_fields[i]->UpdatePhys());
262  m_F[i] = Array< OneD, NekDouble> (m_fields[0]->GetTotPoints(), 0.0);
263  }
264  }
265 
266 
267  /**
268  *
269  */
271  {
272  int nfields = m_fields.num_elements() - 1;
273  for (int k=0 ; k < nfields; ++k)
274  {
275  //Backward Transformation in physical space for time evolution
276  m_fields[k]->BwdTrans_IterPerExp(m_fields[k]->GetCoeffs(),
277  m_fields[k]->UpdatePhys());
278  }
279  }
280 
281  /**
282  *
283  */
285  {
286 
287  int nfields = m_fields.num_elements() - 1;
288  for (int k=0 ; k < nfields; ++k)
289  {
290  //Forward Transformation in physical space for time evolution
291  m_fields[k]->FwdTrans_IterPerExp(m_fields[k]->GetPhys(),m_fields[k]->UpdateCoeffs());
292  }
293  }
294 
295  /**
296  *
297  */
299  {
300  int vVar = m_session->GetVariables().size();
301  Array<OneD, bool> vChecks(vVar, false);
302  vChecks[vVar-1] = true;
303  return vChecks;
304  }
305 
306  /**
307  *
308  */
310  {
311  return m_session->GetVariables().size() - 1;
312  }
313 
314  /**
315  * Explicit part of the method - Advection, Forcing + HOPBCs
316  */
318  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
319  Array<OneD, Array<OneD, NekDouble> > &outarray,
320  const NekDouble time)
321  {
322  EvaluateAdvectionTerms(inarray, outarray);
323 
324  // Smooth advection
326  {
327  for(int i = 0; i < m_nConvectiveFields; ++i)
328  {
329  m_pressure->SmoothField(outarray[i]);
330  }
331  }
332 
333  // Add forcing terms
334  std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
335  for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
336  {
337  (*x)->Apply(m_fields, inarray, outarray, time);
338  }
339 
340  // Calculate High-Order pressure boundary conditions
341  m_extrapolation->EvaluatePressureBCs(inarray,outarray,m_kinvis);
342  }
343 
344  /**
345  * Implicit part of the method - Poisson + nConv*Helmholtz
346  */
348  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
349  Array<OneD, Array<OneD, NekDouble> > &outarray,
350  const NekDouble time,
351  const NekDouble aii_Dt)
352  {
353  // Substep the pressure boundary condition if using substepping
354  m_extrapolation->SubStepSetPressureBCs(inarray,aii_Dt,m_kinvis);
355 
356  // Set up forcing term for pressure Poisson equation
357  SetUpPressureForcing(inarray, m_F, aii_Dt);
358 
359  // Solve Pressure System
360  SolvePressure (m_F[0]);
361 
362  // Set up forcing term for Helmholtz problems
363  SetUpViscousForcing(inarray, m_F, aii_Dt);
364 
365  // Solve velocity system
366  SolveViscous( m_F, outarray, aii_Dt);
367  }
368 
369  /**
370  * Forcing term for Poisson solver solver
371  */
373  const Array<OneD, const Array<OneD, NekDouble> > &fields,
375  const NekDouble aii_Dt)
376  {
377  int i;
378  int physTot = m_fields[0]->GetTotPoints();
379  int nvel = m_velocity.num_elements();
380 
381  m_fields[0]->PhysDeriv(eX,fields[0], Forcing[0]);
382 
383  for(i = 1; i < nvel; ++i)
384  {
385  // Use Forcing[1] as storage since it is not needed for the pressure
386  m_fields[i]->PhysDeriv(DirCartesianMap[i],fields[i],Forcing[1]);
387  Vmath::Vadd(physTot,Forcing[1],1,Forcing[0],1,Forcing[0],1);
388  }
389 
390  Vmath::Smul(physTot,1.0/aii_Dt,Forcing[0],1,Forcing[0],1);
391  }
392 
393  /**
394  * Forcing term for Helmholtz solver
395  */
397  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
399  const NekDouble aii_Dt)
400  {
401  NekDouble aii_dtinv = 1.0/aii_Dt;
402  int phystot = m_fields[0]->GetTotPoints();
403 
404  // Grad p
405  m_pressure->BwdTrans(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
406 
407  int nvel = m_velocity.num_elements();
408  if(nvel == 2)
409  {
410  m_pressure->PhysDeriv(m_pressure->GetPhys(),
411  Forcing[m_velocity[0]],
412  Forcing[m_velocity[1]]);
413  }
414  else
415  {
416  m_pressure->PhysDeriv(m_pressure->GetPhys(),
417  Forcing[m_velocity[0]],
418  Forcing[m_velocity[1]],
419  Forcing[m_velocity[2]]);
420  }
421 
422  // zero convective fields.
423  for(int i = nvel; i < m_nConvectiveFields; ++i)
424  {
425  Vmath::Zero(phystot,Forcing[i],1);
426  }
427 
428  // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
429  // need to be updated for the convected fields.
430  for(int i = 0; i < m_nConvectiveFields; ++i)
431  {
432  Blas::Daxpy(phystot,-aii_dtinv,inarray[i],1,Forcing[i],1);
433  Blas::Dscal(phystot,1.0/m_diffCoeff[i],&(Forcing[i])[0],1);
434  }
435  }
436 
437 
438  /**
439  * Solve pressure system
440  */
443  {
445  // Setup coefficient for equation
446  factors[StdRegions::eFactorLambda] = 0.0;
447 
448  // Solver Pressure Poisson Equation
449  m_pressure->HelmSolve(Forcing, m_pressure->UpdateCoeffs(),
450  NullFlagList, factors);
451 
452  // Add presure to outflow bc if using convective like BCs
453  m_extrapolation->AddPressureToOutflowBCs(m_kinvis);
454  }
455 
456  /**
457  * Solve velocity system
458  */
460  const Array<OneD, const Array<OneD, NekDouble> > &Forcing,
461  Array<OneD, Array<OneD, NekDouble> > &outarray,
462  const NekDouble aii_Dt)
463  {
465 
466  if(m_useSpecVanVisc)
467  {
470  }
471 
472  // Solve Helmholtz system and put in Physical space
473  for(int i = 0; i < m_nConvectiveFields; ++i)
474  {
475  // Setup coefficients for equation
476  factors[StdRegions::eFactorLambda] = 1.0/aii_Dt/m_diffCoeff[i];
477  m_fields[i]->HelmSolve(Forcing[i], m_fields[i]->UpdateCoeffs(),
478  NullFlagList, factors);
479  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),outarray[i]);
480  }
481  }
482 
483 } //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:198
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.
virtual std::string v_GetSubSteppingExtrapolateStr(const std::string &instr)
virtual std::string v_GetExtrapolateStr(void)
void EvaluateAdvectionTerms(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
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:252
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:213
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)
boost::shared_ptr< Equation > EquationSharedPtr
bool m_homogen_dealiasing
Flag to determine if dealiasing is used for homogeneous simulations.
EquationSystemFactory & GetEquationSystemFactory()
Array< OneD, Array< OneD, NekDouble > > m_F
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)
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.
Array< OneD, NekDouble > m_diffCoeff
Diffusion coefficients (will be kinvis for velocities)
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:373
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:299
static FlagList NullFlagList
An empty flag list.
static VarCoeffMap NullVarCoeffMap
Definition: StdRegions.hpp:228
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215