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 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 velcoity fields with homogenous expansion");
117 
118  if(m_useHomo1DSpecVanVisc == false)
119  {
120  m_session->MatchSolverInfo("SpectralVanishingViscosity","True",m_useHomo1DSpecVanVisc,false);
121  }
122 
124  {
125  Array<OneD, unsigned int> planes;
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  {
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 
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 
334  Array<OneD, Array< OneD, NekDouble> > F(m_nConvectiveFields);
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  Vmath::Smul(physTot,1.0/aii_Dt,Forcing[0],1,Forcing[0],1);
393  }
394 
395  /**
396  * Forcing term for Helmholtz solver
397  */
399  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
400  Array<OneD, Array<OneD, NekDouble> > &Forcing,
401  const NekDouble aii_Dt)
402  {
403  NekDouble aii_dtinv = 1.0/aii_Dt;
404  int phystot = m_fields[0]->GetTotPoints();
405 
406  // Grad p
407  m_pressure->BwdTrans(m_pressure->GetCoeffs(),m_pressure->UpdatePhys());
408 
409  int nvel = m_velocity.num_elements();
410  if(nvel == 2)
411  {
412  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1]);
413  }
414  else
415  {
416  m_pressure->PhysDeriv(m_pressure->GetPhys(), Forcing[0], Forcing[1],
417  Forcing[2]);
418  }
419 
420  // Subtract inarray/(aii_dt) and divide by kinvis. Kinvis will
421  // need to be updated for the convected fields.
422  for(int i = 0; i < nvel; ++i)
423  {
424  Blas::Daxpy(phystot,-aii_dtinv,inarray[i],1,Forcing[i],1);
425  Blas::Dscal(phystot,1.0/m_kinvis,&(Forcing[i])[0],1);
426  }
427  }
428 
429 } //end of namespace