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