Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DiffusionLDG.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: DiffusionLDG.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: LDG diffusion class.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <iostream>
38 #include <iomanip>
39 
40 namespace Nektar
41 {
42  namespace SolverUtils
43  {
45  RegisterCreatorFunction("LDG", DiffusionLDG::create);
46 
48  {
49  }
50 
54  {
55  m_session = pSession;
56 
57  m_session->LoadSolverInfo("ShockCaptureType",
58  m_shockCaptureType, "Off");
59 
60  // Setting up the normals
61  int i;
62  int nDim = pFields[0]->GetCoordim(0);
63  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
64 
66  for(i = 0; i < nDim; ++i)
67  {
68  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
69  }
70  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
71  }
72 
74  const int nConvectiveFields,
76  const Array<OneD, Array<OneD, NekDouble> > &inarray,
77  Array<OneD, Array<OneD, NekDouble> > &outarray)
78  {
79  int nBndEdgePts, i, j, k, e;
80  int nDim = fields[0]->GetCoordim(0);
81  int nPts = fields[0]->GetTotPoints();
82  int nCoeffs = fields[0]->GetNcoeffs();
83  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
84 
85  Array<OneD, NekDouble> qcoeffs(nCoeffs);
86  Array<OneD, NekDouble> temp (nCoeffs);
87 
88  Array<OneD, Array<OneD, NekDouble> > fluxvector(nDim);
89 
93 
94  for (j = 0; j < nDim; ++j)
95  {
96  qfield[j] =
97  Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
98  qfieldStd[j] =
99  Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
100  flux[j] =
101  Array<OneD, Array<OneD, NekDouble> >(nConvectiveFields);
102 
103  for (i = 0; i < nConvectiveFields; ++i)
104  {
105  qfield[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
106  qfieldStd[j][i] = Array<OneD, NekDouble>(nPts, 0.0);
107  flux[j][i] = Array<OneD, NekDouble>(nTracePts, 0.0);
108  }
109  }
110 
111  for (k = 0; k < nDim; ++k)
112  {
113  fluxvector[k] = Array<OneD, NekDouble>(nPts, 0.0);
114  }
115 
116  // Compute q_{\eta} and q_{\xi}
117  // Obtain numerical fluxes
118 
119  v_NumFluxforScalar(fields, inarray, flux);
120 
121  for (j = 0; j < nDim; ++j)
122  {
123  for (i = 0; i < nConvectiveFields; ++i)
124  {
125  fields[i]->IProductWRTDerivBase(j, inarray[i], qcoeffs);
126  Vmath::Neg (nCoeffs, qcoeffs, 1);
127  fields[i]->AddTraceIntegral (flux[j][i], qcoeffs);
128  fields[i]->SetPhysState (false);
129  fields[i]->MultiplyByElmtInvMass(qcoeffs, qcoeffs);
130  fields[i]->BwdTrans (qcoeffs, qfield[j][i]);
131  }
132  }
133  // Compute u from q_{\eta} and q_{\xi}
134  // Obtain numerical fluxes
135  v_NumFluxforVector(fields, inarray, qfield, flux[0]);
136 
138  {
139  Array<OneD, NekDouble> muvar(nPts, 0.0);
140  m_ArtificialDiffusionVector(inarray, muvar);
141 
142  int numConvFields = nConvectiveFields;
143 
144  if (m_shockCaptureType == "Smooth")
145  {
146  numConvFields = nConvectiveFields - 1;
147  }
148 
149  for (j = 0; j < nDim; ++j)
150  {
151  for (i = 0; i < numConvFields; ++i)
152  {
153  Vmath::Vmul(nPts,qfield[j][i],1,muvar,1,qfield[j][i],1);
154  }
155  }
156 
157  Array<OneD, NekDouble> FwdMuVar(nTracePts, 0.0);
158  Array<OneD, NekDouble> BwdMuVar(nTracePts, 0.0);
159 
160  fields[0]->GetFwdBwdTracePhys(muvar,FwdMuVar,BwdMuVar);
161 
162  int nBndRegions = fields[0]->GetBndCondExpansions().
163  num_elements();
164  int cnt = 0;
165 
166  for (i = 0; i < nBndRegions; ++i)
167  {
168  // Number of boundary expansion related to that region
169  int nBndEdges = fields[0]->
170  GetBndCondExpansions()[i]->GetExpSize();
171 
172  // Weakly impose boundary conditions by modifying flux
173  // values
174  for (e = 0; e < nBndEdges ; ++e)
175  {
176  nBndEdgePts = fields[0]->GetBndCondExpansions()[i]
177  ->GetExp(e)->GetTotPoints();
178 
179  int id2 = fields[0]->GetTrace()->GetPhys_Offset(
180  fields[0]->GetTraceMap()
181  ->GetBndCondTraceToGlobalTraceMap(cnt++));
182 
183  for (k = 0; k < nBndEdgePts; ++k)
184  {
185  BwdMuVar[id2+k] = 0.0;
186  }
187  }
188  }
189 
190  for(i = 0; i < numConvFields; ++i)
191  {
192  for(k = 0; k < nTracePts; ++k)
193  {
194  flux[0][i][k] =
195  0.5 * (FwdMuVar[k] + BwdMuVar[k]) * flux[0][i][k];
196  }
197  }
198  }
199 
201  Array<OneD, Array<OneD, NekDouble> > qdbase(nDim);
202 
203  for (i = 0; i < nConvectiveFields; ++i)
204  {
205  for (j = 0; j < nDim; ++j)
206  {
207  qdbase[j] = qfield[j][i];
208  }
209  fields[i]->IProductWRTDerivBase(qdbase,tmp);
210 
211  // Evaulate <\phi, \hat{F}\cdot n> - outarray[i]
212  Vmath::Neg (nCoeffs, tmp, 1);
213  fields[i]->AddTraceIntegral (flux[0][i], tmp);
214  fields[i]->SetPhysState (false);
215  fields[i]->MultiplyByElmtInvMass(tmp, tmp);
216  fields[i]->BwdTrans (tmp, outarray[i]);
217  }
218  }
219 
222  const Array<OneD, Array<OneD, NekDouble> > &ufield,
224  {
225  int i, j;
226  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
227  int nvariables = fields.num_elements();
228  int nDim = uflux.num_elements();
229 
230  Array<OneD, NekDouble > Fwd (nTracePts);
231  Array<OneD, NekDouble > Bwd (nTracePts);
232  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
233  Array<OneD, NekDouble > fluxtemp(nTracePts, 0.0);
234 
235  // Get the normal velocity Vn
236  for(i = 0; i < nDim; ++i)
237  {
238  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
239  Vn, 1, Vn, 1);
240  }
241 
242  // Get the sign of (v \cdot n), v = an arbitrary vector
243  // Evaluate upwind flux:
244  // uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
245  for (j = 0; j < nDim; ++j)
246  {
247  for (i = 0; i < nvariables ; ++i)
248  {
249  // Compute Fwd and Bwd value of ufield of i direction
250  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
251 
252  // if Vn >= 0, flux = uFwd, i.e.,
253  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uFwd
254  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uFwd
255 
256  // else if Vn < 0, flux = uBwd, i.e.,
257  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
258  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
259 
260  fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
261  Fwd, Bwd, fluxtemp);
262 
263  // Imposing weak boundary condition with flux
264  // if Vn >= 0, uflux = uBwd at Neumann, i.e.,
265  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick uflux = uBwd
266  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick uflux = uBwd
267 
268  // if Vn >= 0, uflux = uFwd at Neumann, i.e.,
269  // edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uFwd
270  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uFwd
271 
272  if(fields[0]->GetBndCondExpansions().num_elements())
273  {
274  v_WeakPenaltyforScalar(fields, i, ufield[i], fluxtemp);
275  }
276 
277  // if Vn >= 0, flux = uFwd*(tan_{\xi}^- \cdot \vec{n}),
278  // i.e,
279  // edge::eForward, uFwd \(\tan_{\xi}^Fwd \cdot \vec{n})
280  // edge::eBackward, uFwd \(\tan_{\xi}^Bwd \cdot \vec{n})
281 
282  // else if Vn < 0, flux = uBwd*(tan_{\xi}^- \cdot \vec{n}),
283  // i.e,
284  // edge::eForward, uBwd \(\tan_{\xi}^Fwd \cdot \vec{n})
285  // edge::eBackward, uBwd \(\tan_{\xi}^Bwd \cdot \vec{n})
286 
287  Vmath::Vmul(nTracePts,
288  m_traceNormals[j], 1,
289  fluxtemp, 1,
290  uflux[j][i], 1);
291  }
292  }
293  }
294 
295 
296 
299  const int var,
300  const Array<OneD, const NekDouble> &ufield,
301  Array<OneD, NekDouble> &penaltyflux)
302  {
303  int i, e, id1, id2;
304 
305  // Number of boundary regions
306  int nBndEdgePts, nBndEdges;
307  int cnt = 0;
308  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
309  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
310  Array<OneD, NekDouble > uplus(nTracePts);
311 
312  fields[var]->ExtractTracePhys(ufield, uplus);
313  for (i = 0; i < nBndRegions; ++i)
314  {
315  // Number of boundary expansion related to that region
316  nBndEdges = fields[var]->
317  GetBndCondExpansions()[i]->GetExpSize();
318 
319  // Weakly impose boundary conditions by modifying flux values
320  for (e = 0; e < nBndEdges ; ++e)
321  {
322  nBndEdgePts = fields[var]->
323  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
324 
325  id1 = fields[var]->
326  GetBndCondExpansions()[i]->GetPhys_Offset(e);
327 
328  id2 = fields[0]->GetTrace()->
329  GetPhys_Offset(fields[0]->GetTraceMap()->
330  GetBndCondTraceToGlobalTraceMap(cnt++));
331 
332  // For Dirichlet boundary condition: uflux = g_D
333  if (fields[var]->GetBndConditions()[i]->
334  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
335  {
336  Vmath::Vcopy(nBndEdgePts,
337  &(fields[var]->
338  GetBndCondExpansions()[i]->
339  GetPhys())[id1], 1,
340  &penaltyflux[id2], 1);
341  }
342  // For Neumann boundary condition: uflux = u+
343  else if ((fields[var]->GetBndConditions()[i])->
344  GetBoundaryConditionType() == SpatialDomains::eNeumann)
345  {
346  Vmath::Vcopy(nBndEdgePts,
347  &uplus[id2], 1,
348  &penaltyflux[id2], 1);
349  }
350  }
351  }
352  }
353 
354 
355 
358  const Array<OneD, Array<OneD, NekDouble> > &ufield,
361  {
362  int i, j;
363  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
364  int nvariables = fields.num_elements();
365  int nDim = qfield.num_elements();
366 
367  NekDouble C11 = 0.0;
368  Array<OneD, NekDouble > Fwd(nTracePts);
369  Array<OneD, NekDouble > Bwd(nTracePts);
370  Array<OneD, NekDouble > Vn (nTracePts, 0.0);
371  Array<OneD, NekDouble > qFwd (nTracePts);
372  Array<OneD, NekDouble > qBwd (nTracePts);
373  Array<OneD, NekDouble > qfluxtemp(nTracePts, 0.0);
374  Array<OneD, NekDouble > uterm(nTracePts);
375  /*
376  // Setting up the normals
377  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDim);
378  for(i = 0; i < nDim; ++i)
379  {
380  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
381  }
382  fields[0]->GetTrace()->GetNormals(m_traceNormals);
383  */
384 
385  // Get the normal velocity Vn
386  for(i = 0; i < nDim; ++i)
387  {
388  Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
389  Vn, 1, Vn, 1);
390  }
391  // Evaulate upwind flux:
392  // qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
393  for (i = 0; i < nvariables; ++i)
394  {
395  qflux[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
396  for (j = 0; j < nDim; ++j)
397  {
398  // Compute Fwd and Bwd value of ufield of jth direction
399  fields[i]->GetFwdBwdTracePhys(qfield[j][i], qFwd, qBwd);
400 
401  // if Vn >= 0, flux = uFwd, i.e.,
402  // edge::eForward, if V*n>=0 <=> V*n_F>=0, pick
403  // qflux = qBwd = q+
404  // edge::eBackward, if V*n>=0 <=> V*n_B<0, pick
405  // qflux = qBwd = q-
406 
407  // else if Vn < 0, flux = uBwd, i.e.,
408  // edge::eForward, if V*n<0 <=> V*n_F<0, pick
409  // qflux = qFwd = q-
410  // edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
411  // qflux = qFwd = q+
412 
413  fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
414  qBwd, qFwd,
415  qfluxtemp);
416 
417  Vmath::Vmul(nTracePts,
418  m_traceNormals[j], 1,
419  qfluxtemp, 1,
420  qfluxtemp, 1);
421 
422  // Generate Stability term = - C11 ( u- - u+ )
423  fields[i]->GetFwdBwdTracePhys(ufield[i], Fwd, Bwd);
424 
425  Vmath::Vsub(nTracePts,
426  Fwd, 1, Bwd, 1,
427  uterm, 1);
428 
429  Vmath::Smul(nTracePts,
430  -1.0 * C11, uterm, 1,
431  uterm, 1);
432 
433  // Flux = {Fwd, Bwd} * (nx, ny, nz) + uterm * (nx, ny)
434  Vmath::Vadd(nTracePts,
435  uterm, 1,
436  qfluxtemp, 1,
437  qfluxtemp, 1);
438 
439  // Imposing weak boundary condition with flux
440  if (fields[0]->GetBndCondExpansions().num_elements())
441  {
442  v_WeakPenaltyforVector(fields, i, j,
443  qfield[j][i],
444  qfluxtemp, C11);
445  }
446 
447  // q_hat \cdot n = (q_xi \cdot n_xi) or (q_eta \cdot n_eta)
448  // n_xi = n_x * tan_xi_x + n_y * tan_xi_y + n_z * tan_xi_z
449  // n_xi = n_x * tan_eta_x + n_y * tan_eta_y + n_z*tan_eta_z
450  Vmath::Vadd(nTracePts,
451  qfluxtemp, 1,
452  qflux[i], 1,
453  qflux[i], 1);
454  }
455  }
456  }
457 
458 
459 
460  /**
461  * Diffusion: Imposing weak boundary condition for q with flux
462  * uflux = g_D on Dirichlet boundary condition
463  * uflux = u_Fwd on Neumann boundary condition
464  */
467  const int var,
468  const int dir,
469  const Array<OneD, const NekDouble> &qfield,
470  Array<OneD, NekDouble> &penaltyflux,
471  NekDouble C11)
472  {
473  int i, e, id1, id2;
474  int nBndEdges, nBndEdgePts;
475  int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
476  int nTracePts = fields[0]->GetTrace()->GetTotPoints();
477 
478  Array<OneD, NekDouble > uterm(nTracePts);
479  Array<OneD, NekDouble > qtemp(nTracePts);
480  int cnt = 0;
481 
482  /*
483  // Setting up the normals
484  m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(nDim);
485  for(i = 0; i < nDim; ++i)
486  {
487  m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
488  }
489  fields[0]->GetTrace()->GetNormals(m_traceNormals);
490  */
491 
492  fields[var]->ExtractTracePhys(qfield, qtemp);
493 
494  for (i = 0; i < nBndRegions; ++i)
495  {
496  nBndEdges = fields[var]->
497  GetBndCondExpansions()[i]->GetExpSize();
498 
499  // Weakly impose boundary conditions by modifying flux values
500  for (e = 0; e < nBndEdges ; ++e)
501  {
502  nBndEdgePts = fields[var]->
503  GetBndCondExpansions()[i]->GetExp(e)->GetTotPoints();
504 
505  id1 = fields[var]->
506  GetBndCondExpansions()[i]->GetPhys_Offset(e);
507 
508  id2 = fields[0]->GetTrace()->
509  GetPhys_Offset(fields[0]->GetTraceMap()->
510  GetBndCondTraceToGlobalTraceMap(cnt++));
511 
512  // For Dirichlet boundary condition:
513  //qflux = q+ - C_11 (u+ - g_D) (nx, ny)
514  if(fields[var]->GetBndConditions()[i]->
515  GetBoundaryConditionType() == SpatialDomains::eDirichlet)
516  {
517  Vmath::Vmul(nBndEdgePts,
518  &m_traceNormals[dir][id2], 1,
519  &qtemp[id2], 1,
520  &penaltyflux[id2], 1);
521  }
522  // For Neumann boundary condition: qflux = g_N
523  else if((fields[var]->GetBndConditions()[i])->
524  GetBoundaryConditionType() == SpatialDomains::eNeumann)
525  {
526  Vmath::Vmul(nBndEdgePts,
527  &m_traceNormals[dir][id2], 1,
528  &(fields[var]->
529  GetBndCondExpansions()[i]->
530  GetPhys())[id1], 1,
531  &penaltyflux[id2], 1);
532  }
533  }
534  }
535  }
536 
537  }
538 }
DiffusionArtificialDiffusion m_ArtificialDiffusionVector
Definition: Diffusion.h:135
DiffusionFactory & GetDiffusionFactory()
Definition: Diffusion.cpp:42
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:471
static DiffusionSharedPtr create(std::string diffType)
Definition: DiffusionLDG.h:48
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
virtual void v_NumFluxforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &uflux)
LibUtilities::SessionReaderSharedPtr m_session
Definition: DiffusionLDG.h:61
virtual void v_WeakPenaltyforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const int dir, const Array< OneD, const NekDouble > &qfield, Array< OneD, NekDouble > &penaltyflux, NekDouble C11)
virtual void v_NumFluxforVector(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &ufield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &qfield, Array< OneD, Array< OneD, NekDouble > > &qflux)
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
virtual void v_Diffuse(const int nConvective, const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:329
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Definition: DiffusionLDG.h:60
virtual void v_WeakPenaltyforScalar(const Array< OneD, MultiRegions::ExpListSharedPtr > &fields, const int var, const Array< OneD, const NekDouble > &ufield, Array< OneD, NekDouble > &penaltyflux)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
virtual void v_InitObject(LibUtilities::SessionReaderSharedPtr pSession, Array< OneD, MultiRegions::ExpListSharedPtr > pFields)
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
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:169