Nektar++
MMFDiffusion.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File MMFDiffusion.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: MMFDiffusion.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <iomanip>
36 #include <iostream>
37 
38 #include <boost/algorithm/string.hpp>
39 #include <boost/core/ignore_unused.hpp>
40 
41 #include <DiffusionSolver/EquationSystems/MMFDiffusion.h>
42 #include <LibUtilities/BasicUtils/SessionReader.h>
43 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
44 #include <SolverUtils/Driver.h>
45 
46 #include <boost/math/special_functions/spherical_harmonic.hpp>
47 using namespace std;
48 using namespace Nektar::SolverUtils;
49 using namespace Nektar;
50 
51 namespace Nektar
52 {
53 string MMFDiffusion::className =
54  SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
55  "MMFDiffusion", MMFDiffusion::create, "MMFDiffusion equation.");
56 
57 MMFDiffusion::MMFDiffusion(const LibUtilities::SessionReaderSharedPtr &pSession,
58  const SpatialDomains::MeshGraphSharedPtr &pGraph)
59  : UnsteadySystem(pSession, pGraph), MMFSystem(pSession, pGraph)
60 {
61 }
62 
63 void MMFDiffusion::v_InitObject(bool DeclareFields)
64 {
65  UnsteadySystem::v_InitObject(DeclareFields);
66 
67  int nq = m_fields[0]->GetNpoints();
68  int nvar = m_fields.size();
69  int MFdim = 3;
70 
71  // Diffusivity coefficient for e^j
72  m_epsilon = Array<OneD, NekDouble>(MFdim);
73  m_session->LoadParameter("epsilon0", m_epsilon[0], 1.0);
74  m_session->LoadParameter("epsilon1", m_epsilon[1], 1.0);
75  m_session->LoadParameter("epsilon2", m_epsilon[2], 1.0);
76 
77  // Diffusivity coefficient for u^j
78  m_epsu = Array<OneD, NekDouble>(nvar + 1);
79  m_session->LoadParameter("epsu0", m_epsu[0], 1.0);
80  m_session->LoadParameter("epsu1", m_epsu[1], 1.0);
81 
82  m_session->LoadParameter("InitPtx", m_InitPtx, 0.0);
83  m_session->LoadParameter("InitPty", m_InitPty, 0.0);
84  m_session->LoadParameter("InitPtz", m_InitPtz, 0.0);
85 
86  int shapedim = m_fields[0]->GetShapeDimension();
87  Array<OneD, Array<OneD, NekDouble>> Anisotropy(shapedim);
88  for (int j = 0; j < shapedim; ++j)
89  {
90  Anisotropy[j] = Array<OneD, NekDouble>(nq, 1.0);
91  Vmath::Fill(nq, sqrt(m_epsilon[j]), &Anisotropy[j][0], 1);
92  }
93 
94  MMFSystem::MMFInitObject(Anisotropy);
95 
96  // Define ProblemType
97  if (m_session->DefinesSolverInfo("TESTTYPE"))
98  {
99  std::string TestTypeStr = m_session->GetSolverInfo("TESTTYPE");
100  int i;
101  for (i = 0; i < (int)SIZE_TestType; ++i)
102  {
103  if (boost::iequals(TestTypeMap[i], TestTypeStr))
104  {
105  m_TestType = (TestType)i;
106  break;
107  }
108  }
109  }
110  else
111  {
112  m_TestType = (TestType)0;
113  }
114 
115  if (m_session->DefinesSolverInfo("INITWAVETYPE"))
116  {
117  std::string InitWaveTypeStr = m_session->GetSolverInfo("INITWAVETYPE");
118  for (int i = 0; i < (int)SIZE_TestType; ++i)
119  {
120  if (boost::iequals(InitWaveTypeMap[i], InitWaveTypeStr))
121  {
122  m_InitWaveType = (InitWaveType)i;
123  break;
124  }
125  }
126  }
127  else
128  {
129  m_InitWaveType = (InitWaveType)0;
130  }
131 
132  StdRegions::VarCoeffType MMFCoeffs[15] = {
133  StdRegions::eVarCoeffMF1x, StdRegions::eVarCoeffMF1y,
134  StdRegions::eVarCoeffMF1z, StdRegions::eVarCoeffMF1Div,
135  StdRegions::eVarCoeffMF1Mag, StdRegions::eVarCoeffMF2x,
136  StdRegions::eVarCoeffMF2y, StdRegions::eVarCoeffMF2z,
137  StdRegions::eVarCoeffMF2Div, StdRegions::eVarCoeffMF2Mag,
138  StdRegions::eVarCoeffMF3x, StdRegions::eVarCoeffMF3y,
139  StdRegions::eVarCoeffMF3z, StdRegions::eVarCoeffMF3Div,
140  StdRegions::eVarCoeffMF3Mag};
141 
142  int indx;
143  Array<OneD, NekDouble> tmp(nq);
144  for (int k = 0; k < MFdim; ++k)
145  {
146  // For Moving Frames
147  indx = 5 * k;
148 
149  for (int j = 0; j < m_spacedim; ++j)
150  {
151  m_varcoeff[MMFCoeffs[indx + j]] = Array<OneD, NekDouble>(nq, 0.0);
152  Vmath::Vcopy(nq, &m_movingframes[k][j * nq], 1,
153  &m_varcoeff[MMFCoeffs[indx + j]][0], 1);
154  }
155 
156  // m_DivMF
157  m_varcoeff[MMFCoeffs[indx + 3]] = Array<OneD, NekDouble>(nq, 0.0);
158  Vmath::Vcopy(nq, &m_DivMF[k][0], 1, &m_varcoeff[MMFCoeffs[indx + 3]][0],
159  1);
160 
161  // \| e^k \|
162  m_varcoeff[MMFCoeffs[indx + 4]] = Array<OneD, NekDouble>(nq, 0.0);
163  tmp = Array<OneD, NekDouble>(nq, 0.0);
164  for (int i = 0; i < m_spacedim; ++i)
165  {
166  Vmath::Vvtvp(nq, &m_movingframes[k][i * nq], 1,
167  &m_movingframes[k][i * nq], 1, &tmp[0], 1, &tmp[0], 1);
168  }
169 
170  Vmath::Vcopy(nq, &tmp[0], 1, &m_varcoeff[MMFCoeffs[indx + 4]][0], 1);
171  }
172 
173  if (!m_explicitDiffusion)
174  {
175  m_ode.DefineImplicitSolve(&MMFDiffusion::DoImplicitSolve, this);
176  }
177 
178  m_ode.DefineOdeRhs(&MMFDiffusion::DoOdeRhs, this);
179 }
180 
181 /**
182  *
183  */
184 MMFDiffusion::~MMFDiffusion()
185 {
186 }
187 
188 /**OdeRhs
189  * @param inarray Input array.
190  * @param outarray Output array.
191  * @param time Current simulation time.
192  * @param lambda Timestep.
193  */
194 void MMFDiffusion::DoImplicitSolve(
195  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
196  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time,
197  const NekDouble lambda)
198 {
199  int nvariables = inarray.size();
200  int nq = m_fields[0]->GetNpoints();
201 
202  StdRegions::ConstFactorMap factors;
203  factors[StdRegions::eFactorTau] = 1.0;
204 
205  Array<OneD, Array<OneD, NekDouble>> F(nvariables);
206  factors[StdRegions::eFactorLambda] = 1.0 / lambda;
207  F[0] = Array<OneD, NekDouble>(nq * nvariables);
208 
209  for (int n = 1; n < nvariables; ++n)
210  {
211  F[n] = F[n - 1] + nq;
212  // cout << "F["<< n<<"=" << F[n][1] <<endl;
213  }
214 
215  // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
216  // inarray = input: \hat{rhs} -> output: \hat{Y}
217  // outarray = output: nabla^2 \hat{Y}
218  // where \hat = modal coeffs
219  SetBoundaryConditions(time);
220 
221  for (int i = 0; i < nvariables; ++i)
222  {
223  factors[StdRegions::eFactorLambda] = 1.0 / lambda / m_epsu[i];
224 
225  // Multiply 1.0/timestep
226  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda], inarray[i], 1,
227  F[i], 1);
228 
229  /* for (int k = 0; k < 15; ++k)
230  cout << "inarray["<<i << "]"<< k<<"=" << inarray[i][k]<<endl;*/
231  // Solve a system of equations with Helmholtz solver and transform
232  // back into physical space.
233  m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(), factors,
234  m_varcoeff);
235 
236  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
237  /* Array<OneD, NekDouble> coefarray = m_fields[i]->GetCoeffs();
238  for (int k = 0; k < 15; ++k)
239  cout << "inarray["<< k<<"=" << coefarray[k]<<endl;*/
240  }
241  /* for (int kk = 0; kk < 15; ++kk)
242  cout << "inarray["<< kk<<"=" <<
243  m_varcoeff[StdRegions::eVarCoeffMF3Mag][kk]<<endl;*/
244 }
245 
246 /**
247  *
248  */
249 void MMFDiffusion::DoOdeRhs(
250  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
251  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
252 {
253  int nq = GetTotPoints();
254 
255  switch (m_TestType)
256  {
257  case eTestPlane:
258  {
259 
260  Array<OneD, NekDouble> x(nq);
261  Array<OneD, NekDouble> y(nq);
262  Array<OneD, NekDouble> z(nq);
263 
264  m_fields[0]->GetCoords(x, y, z);
265 
266  for (int k = 0; k < nq; k++)
267  {
268  outarray[0][k] = (m_epsilon[0] + m_epsilon[1] - 1.0) * m_pi *
269  m_pi * exp(-1.0 * m_pi * m_pi * time) *
270  sin(m_pi * x[k]) * cos(m_pi * y[k]);
271  }
272  }
273  break;
274 
275  case eTestCube:
276  {
277 
278  Array<OneD, NekDouble> x(nq);
279  Array<OneD, NekDouble> y(nq);
280  Array<OneD, NekDouble> z(nq);
281 
282  m_fields[0]->GetCoords(x, y, z);
283 
284  for (int k = 0; k < nq; k++)
285  {
286  outarray[0][k] =
287  (m_epsilon[0] + m_epsilon[1] + m_epsilon[2] - 1.0) * m_pi *
288  m_pi * exp(-1.0 * m_pi * m_pi * time) * sin(m_pi * x[k]) *
289  sin(m_pi * y[k]) * sin(m_pi * z[k]);
290  }
291  }
292  break;
293 
294  case eTestLinearSphere:
295  {
296  Array<OneD, NekDouble> temp(nq);
297 
298  NekDouble A = 2.0;
299  NekDouble B = 5.0;
300 
301  NekDouble m_a, m_b, m_c, m_d;
302  m_a = B - 1.0;
303  m_b = A * A;
304  m_c = -1.0 * B;
305  m_d = -1.0 * A * A;
306 
307  temp = Array<OneD, NekDouble>(nq, 0.0);
308  Vmath::Svtvp(nq, m_a, &inarray[0][0], 1, &temp[0], 1, &temp[0], 1);
309  Vmath::Svtvp(nq, m_b, &inarray[1][0], 1, &temp[0], 1,
310  &outarray[0][0], 1);
311 
312  temp = Array<OneD, NekDouble>(nq, 0.0);
313  Vmath::Svtvp(nq, m_c, &inarray[0][0], 1, &temp[0], 1, &temp[0], 1);
314  Vmath::Svtvp(nq, m_d, &inarray[1][0], 1, &temp[0], 1,
315  &outarray[1][0], 1);
316  }
317  break;
318 
319  case eTestNonlinearSphere:
320  {
321  NekDouble A = 2.0;
322  NekDouble B = 5.0;
323 
324  Array<OneD, NekDouble> Aonevec(nq, A);
325 
326  // cube = phys0*phys0*phy1
327  Array<OneD, NekDouble> cube(nq);
328  Vmath::Vmul(nq, &inarray[0][0], 1, &inarray[0][0], 1, &cube[0], 1);
329  Vmath::Vmul(nq, &inarray[1][0], 1, &cube[0], 1, &cube[0], 1);
330 
331  // outarray[0] = A - B*phy0 + phy0*phy0*phy1 - phy0
332  NekDouble coeff = -1.0 * B - 1.0;
333  Array<OneD, NekDouble> tmp(nq);
334  Vmath::Svtvp(nq, coeff, &inarray[0][0], 1, &cube[0], 1, &tmp[0], 1);
335  Vmath::Vadd(nq, &Aonevec[0], 1, &tmp[0], 1, &outarray[0][0], 1);
336 
337  // outarray[1] = B*phys0 - phy0*phy0*phy1
338  Vmath::Svtvm(nq, B, &inarray[0][0], 1, &cube[0], 1, &outarray[1][0],
339  1);
340  }
341  break;
342 
343  case eFHNStandard:
344  {
345  // \phi - \phi^3/3 - \psi
346  NekDouble a = 0.12;
347  NekDouble b = 0.011;
348  NekDouble c1 = 0.175;
349  NekDouble c2 = 0.03;
350  NekDouble d = 0.55;
351 
352  Array<OneD, NekDouble> tmp(nq);
353 
354  // Reaction for \phi = c1 \phi ( \phi - a)*(1 - \phi) - c2 v
355  Vmath::Smul(nq, -1.0 * c1, inarray[0], 1, outarray[0], 1);
356  Vmath::Sadd(nq, -1.0 * a, inarray[0], 1, tmp, 1);
357  Vmath::Vmul(nq, tmp, 1, inarray[0], 1, outarray[0], 1);
358  Vmath::Sadd(nq, -1.0, inarray[0], 1, tmp, 1);
359  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
360 
361  Vmath::Smul(nq, -1.0 * c2, inarray[1], 1, tmp, 1);
362  Vmath::Vadd(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
363 
364  // Reaction for \psi = b (\phi - d \psi )
365  Vmath::Svtvp(nq, -1.0 * d, inarray[1], 1, inarray[0], 1,
366  outarray[1], 1);
367  Vmath::Smul(nq, b, outarray[1], 1, outarray[1], 1);
368  }
369  break;
370 
371  case eFHNRogers:
372  {
373  NekDouble a = 0.13;
374  NekDouble b = 0.013;
375  NekDouble c1 = 0.26;
376  NekDouble c2 = 0.1;
377  NekDouble d = 1.0;
378 
379  Array<OneD, NekDouble> tmp(nq);
380 
381  // Reaction for \phi = c1 \phi ( \phi - a)*(1 - \phi) - c2 u v
382  Vmath::Smul(nq, -1.0 * c1, inarray[0], 1, outarray[0], 1);
383  Vmath::Sadd(nq, -1.0 * a, inarray[0], 1, tmp, 1);
384  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
385  Vmath::Sadd(nq, -1.0, inarray[0], 1, tmp, 1);
386  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
387 
388  Vmath::Vmul(nq, inarray[0], 1, inarray[1], 1, tmp, 1);
389  Vmath::Smul(nq, -1.0 * c2, tmp, 1, tmp, 1);
390  Vmath::Vadd(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
391 
392  // Reaction for \psi = b (\phi - d \psi )
393  Vmath::Svtvp(nq, -1.0 * d, inarray[1], 1, inarray[0], 1,
394  outarray[1], 1);
395  Vmath::Smul(nq, b, outarray[1], 1, outarray[1], 1);
396  }
397  break;
398 
399  case eFHNAlievPanf:
400  {
401 
402  NekDouble a = 0.15;
403  NekDouble c1 = 8.0;
404  NekDouble c2 = 1.0;
405  NekDouble c0 = 0.002;
406  NekDouble mu1 = 0.2;
407  NekDouble mu2 = 0.3;
408 
409  Array<OneD, NekDouble> tmp(nq);
410 
411  // Reaction for \phi = c1 \phi ( \phi - a)*(1 - \phi) - c2 u v
412  Vmath::Smul(nq, -1.0 * c1, inarray[0], 1, outarray[0], 1);
413  Vmath::Sadd(nq, -1.0 * a, inarray[0], 1, tmp, 1);
414  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
415  Vmath::Sadd(nq, -1.0, inarray[0], 1, tmp, 1);
416  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
417 
418  Vmath::Vmul(nq, inarray[0], 1, inarray[1], 1, tmp, 1);
419  Vmath::Smul(nq, -1.0 * c2, tmp, 1, tmp, 1);
420  Vmath::Vadd(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
421 
422  // Reaction for \psi = (c0 + (\mu1 \psi/(\mu2+\phi) ) )*(-\psi - c1
423  // * \phi*(\phi - a - 1) )
424 
425  Vmath::Smul(nq, mu1, inarray[1], 1, outarray[1], 1);
426  Vmath::Sadd(nq, mu2, inarray[0], 1, tmp, 1);
427  Vmath::Vdiv(nq, outarray[1], 1, tmp, 1, outarray[1], 1);
428  Vmath::Sadd(nq, c0, outarray[1], 1, outarray[1], 1);
429 
430  Vmath::Sadd(nq, (-a - 1.0), inarray[0], 1, tmp, 1);
431  Vmath::Vmul(nq, inarray[0], 1, tmp, 1, tmp, 1);
432  Vmath::Smul(nq, c1, tmp, 1, tmp, 1);
433  Vmath::Vadd(nq, inarray[1], 1, tmp, 1, tmp, 1);
434  Vmath::Neg(nq, tmp, 1);
435 
436  Vmath::Vmul(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
437  }
438  break;
439 
440  default:
441  break;
442  }
443 }
444 
445 /**
446  *
447  */
448 void MMFDiffusion::v_SetInitialConditions(NekDouble initialtime,
449  bool dumpInitialConditions,
450  const int domain)
451 {
452  boost::ignore_unused(domain);
453 
454  int nq = GetTotPoints();
455 
456  switch (m_TestType)
457  {
458  case eTestPlane:
459  {
460  Array<OneD, NekDouble> u(nq);
461 
462  TestPlaneProblem(initialtime, u);
463  m_fields[0]->SetPhys(u);
464  }
465  break;
466 
467  case eTestCube:
468  {
469  Array<OneD, NekDouble> u(nq);
470 
471  TestCubeProblem(initialtime, u);
472  m_fields[0]->SetPhys(u);
473  /*for (int k=0; k<nq; ++k)
474  {
475  //for (int j=0; j<m_spacedim; ++j)
476  //{
477  cout << "_varcoeff" << u[k] <<endl;
478  // }
479  }*/
480  }
481  break;
482 
483  case eTestLinearSphere:
484  case eTestNonlinearSphere:
485  {
486  Array<OneD, NekDouble> u(nq);
487  Array<OneD, NekDouble> v(nq);
488 
489  Morphogenesis(initialtime, 0, u);
490  Morphogenesis(initialtime, 1, v);
491 
492  m_fields[0]->SetPhys(u);
493  m_fields[1]->SetPhys(v);
494  }
495  break;
496 
497  case eFHNStandard:
498  case eFHNRogers:
499  case eFHNAlievPanf:
500  {
501  Array<OneD, NekDouble> Zero(nq, 0.0);
502  m_fields[0]->SetPhys(PlanePhiWave());
503  m_fields[1]->SetPhys(Zero);
504  }
505  break;
506 
507  default:
508  {
509  EquationSystem::v_SetInitialConditions(initialtime, false);
510  }
511  break;
512  }
513 
514  // forward transform to fill the modal coeffs
515  for (int i = 0; i < m_fields.size(); ++i)
516  {
517  m_fields[i]->SetPhysState(true);
518  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
519  m_fields[i]->UpdateCoeffs());
520  }
521 
522  if (dumpInitialConditions)
523  {
524  std::string outname = m_sessionName + "_initial.chk";
525  WriteFld(outname);
526  }
527 }
528 
529 void MMFDiffusion::TestPlaneProblem(const NekDouble time,
530  Array<OneD, NekDouble> &outfield)
531 
532 {
533  int nq = GetTotPoints();
534 
535  Array<OneD, NekDouble> x(nq);
536  Array<OneD, NekDouble> y(nq);
537  Array<OneD, NekDouble> z(nq);
538 
539  m_fields[0]->GetCoords(x, y, z);
540 
541  outfield = Array<OneD, NekDouble>(nq);
542  for (int k = 0; k < nq; k++)
543  {
544  outfield[k] = exp(-1.0 * m_pi * m_pi * time) * sin(m_pi * x[k]) *
545  cos(m_pi * y[k]);
546  }
547 }
548 
549 void MMFDiffusion::TestCubeProblem(const NekDouble time,
550  Array<OneD, NekDouble> &outfield)
551 
552 {
553  int nq = GetTotPoints();
554 
555  Array<OneD, NekDouble> x(nq);
556  Array<OneD, NekDouble> y(nq);
557  Array<OneD, NekDouble> z(nq);
558 
559  m_fields[0]->GetCoords(x, y, z);
560 
561  outfield = Array<OneD, NekDouble>(nq);
562  for (int k = 0; k < nq; k++)
563  {
564  outfield[k] = exp(-1.0 * m_pi * m_pi * time) * sin(m_pi * x[k]) *
565  sin(m_pi * y[k]) * sin(m_pi * z[k]);
566  }
567 }
568 
569 void MMFDiffusion::Morphogenesis(const NekDouble time, unsigned int field,
570  Array<OneD, NekDouble> &outfield)
571 {
572  int nq = GetTotPoints();
573 
574  int i, m, n, ind;
575  NekDouble a_n, d_n, gamma_n;
576  NekDouble A_mn, C_mn, theta, phi, radius;
577 
578  std::complex<double> Spericharmonic, delta_n, temp;
579  std::complex<double> varphi0, varphi1;
580  std::complex<double> B_mn, D_mn;
581 
582  // Set some parameter values
583  int Maxn = 6;
584  int Maxm = 2 * Maxn - 1;
585 
586  NekDouble A = 2.0;
587  NekDouble B = 5.0;
588 
589  NekDouble m_mu = 0.001;
590  NekDouble m_nu = 0.002;
591 
592  NekDouble m_a, m_b, m_c, m_d;
593 
594  m_a = B - 1.0;
595  m_b = A * A;
596  m_c = -1.0 * B;
597  m_d = -1.0 * A * A;
598 
599  Array<OneD, Array<OneD, NekDouble>> Ainit(Maxn);
600  Array<OneD, Array<OneD, NekDouble>> Binit(Maxn);
601 
602  for (i = 0; i < Maxn; ++i)
603  {
604  Ainit[i] = Array<OneD, NekDouble>(Maxm, 0.0);
605  Binit[i] = Array<OneD, NekDouble>(Maxm, 0.0);
606  }
607 
608  Ainit[5][0] = -0.5839;
609  Ainit[5][1] = -0.8436;
610  Ainit[5][2] = -0.4764;
611  Ainit[5][3] = 0.6475;
612  Ainit[5][4] = 0.1886;
613  Ainit[5][5] = 0.8709;
614  Ainit[5][6] = -0.8338;
615  Ainit[5][7] = 0.1795;
616  Ainit[5][8] = -0.7873;
617  Ainit[5][9] = 0.8842;
618  Ainit[5][10] = 0.2943;
619 
620  Binit[5][0] = -0.6263;
621  Binit[5][1] = 0.9803;
622  Binit[5][2] = 0.7222;
623  Binit[5][3] = 0.5945;
624  Binit[5][4] = 0.6026;
625  Binit[5][5] = -0.2076;
626  Binit[5][6] = 0.4556;
627  Binit[5][7] = 0.6024;
628  Binit[5][8] = 0.9695;
629  Binit[5][9] = -0.4936;
630  Binit[5][10] = 0.1098;
631 
632  Array<OneD, NekDouble> u(nq);
633  Array<OneD, NekDouble> v(nq);
634  Array<OneD, NekDouble> x(nq);
635  Array<OneD, NekDouble> y(nq);
636  Array<OneD, NekDouble> z(nq);
637 
638  m_fields[0]->GetCoords(x, y, z);
639  for (int i = 0; i < nq; ++i)
640  {
641  radius = sqrt(x[i] * x[i] + y[i] * y[i] + z[i] * z[i]);
642 
643  // theta is in [0, pi]
644  theta = asin(z[i] / radius) + 0.5 * m_pi;
645 
646  // phi is in [0, 2*pi]
647  phi = atan2(y[i], x[i]) + m_pi;
648 
649  varphi0 = 0.0 * varphi0;
650  varphi1 = 0.0 * varphi1;
651  for (n = 0; n < Maxn; ++n)
652  {
653  // Set up parameters
654  a_n = m_a - m_mu * (n * (n + 1) / radius / radius);
655  d_n = m_d - m_nu * (n * (n + 1) / radius / radius);
656 
657  gamma_n = 0.5 * (a_n + d_n);
658 
659  temp = (a_n + d_n) * (a_n + d_n) - 4.0 * (a_n * d_n - m_b * m_c);
660  delta_n = 0.5 * sqrt(temp);
661 
662  for (m = -n; m <= n; ++m)
663  {
664  ind = m + n;
665  A_mn = Ainit[n][ind];
666  C_mn = Binit[n][ind];
667 
668  B_mn = ((a_n - gamma_n) * Ainit[n][ind] + m_b * Binit[n][ind]) /
669  delta_n;
670  D_mn = (m_c * Ainit[n][ind] + (d_n - gamma_n) * Binit[n][ind]) /
671  delta_n;
672 
673  Spericharmonic =
674  boost::math::spherical_harmonic(n, m, theta, phi);
675  varphi0 += exp(gamma_n * time) *
676  (A_mn * cosh(delta_n * time) +
677  B_mn * sinh(delta_n * time)) *
678  Spericharmonic;
679  varphi1 += exp(gamma_n * time) *
680  (C_mn * cosh(delta_n * time) +
681  D_mn * sinh(delta_n * time)) *
682  Spericharmonic;
683  }
684  }
685 
686  u[i] = varphi0.real();
687  v[i] = varphi1.real();
688  }
689 
690  switch (field)
691  {
692  case 0:
693  {
694  outfield = u;
695  }
696  break;
697 
698  case 1:
699  {
700  outfield = v;
701  }
702  break;
703  }
704 }
705 
706 Array<OneD, NekDouble> MMFDiffusion::PlanePhiWave()
707 {
708  int nq = GetTotPoints();
709  Array<OneD, NekDouble> outarray(nq, 0.0);
710 
711  Array<OneD, NekDouble> x(nq);
712  Array<OneD, NekDouble> y(nq);
713  Array<OneD, NekDouble> z(nq);
714 
715  m_fields[0]->GetCoords(x, y, z);
716 
717  NekDouble xmin, ymin, xmax;
718 
719  xmin = Vmath::Vmin(nq, x, 1);
720  xmax = Vmath::Vmax(nq, x, 1);
721  ymin = Vmath::Vmin(nq, y, 1);
722 
723  NekDouble xp, yp, xp2;
724  for (int i = 0; i < nq; i++)
725  {
726  switch (m_InitWaveType)
727  {
728  case eLeft:
729  {
730  NekDouble radiusofinit = 4.0;
731  NekDouble frontstiff = 0.1;
732 
733  xp = x[i] - xmin;
734  outarray[i] =
735  1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff));
736  }
737  break;
738 
739  case eBothEnds:
740  {
741  NekDouble radiusofinit = 3.0;
742  NekDouble frontstiff = 0.1;
743 
744  xp = x[i] - xmin;
745  xp2 = x[i] - xmax;
746 
747  outarray[i] =
748  1.0 / (1.0 +
749  exp((sqrt(xp * xp) - radiusofinit) / frontstiff)) +
750  1.0 / (1.0 +
751  exp((sqrt(xp2 * xp2) - radiusofinit) / frontstiff));
752  }
753  break;
754 
755  case eCenter:
756  {
757  NekDouble radiusofinit = 6.0;
758  NekDouble frontstiff = 0.1;
759 
760  // NekDouble xc = 0.5*(Vmath::Vmax(nq, x, 1) + Vmath::Vmin(nq,
761  // x, 1));
762 
763  xp = x[i] - xmin;
764  outarray[i] =
765  1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff));
766  }
767  break;
768 
769  case eLeftBottomCorner:
770  {
771  NekDouble radiusofinit = 6.0;
772  NekDouble frontstiff = 0.1;
773  NekDouble bs = 2.0;
774 
775  xp = x[i] - xmin;
776  yp = y[i] - ymin;
777  outarray[i] =
778  1.0 /
779  (1.0 + exp((sqrt(xp * xp + yp * yp) / bs - radiusofinit) /
780  frontstiff));
781  }
782  break;
783 
784  case ePoint:
785  {
786  NekDouble xloc, yloc, zloc, rad;
787  NekDouble radiusofinit = 10.0;
788 
789  xloc = x[i] - m_InitPtx;
790  yloc = y[i] - m_InitPty;
791  zloc = z[i] - m_InitPtz;
792 
793  rad = sqrt(xloc * xloc + yloc * yloc + zloc * zloc);
794 
795  xloc = xloc / radiusofinit;
796  yloc = yloc / radiusofinit;
797  zloc = zloc / radiusofinit;
798 
799  if (rad < radiusofinit)
800  {
801  outarray[i] =
802  exp(-(1.0 / 2.0) *
803  (xloc * xloc + yloc * yloc + zloc * zloc));
804  }
805 
806  else
807  {
808  outarray[i] = 0.0;
809  }
810  }
811  break;
812 
813  case eSpiralDock:
814  {
815  NekDouble radiusofinit = 3.0;
816  NekDouble frontstiff = 0.1;
817  xp = x[i] - 4.0;
818  yp = y[i];
819  outarray[i] =
820  (1.0 / (1.0 + exp(2.0 * yp))) *
821  (1.0 / (1.0 + exp(-2.0 * xp))) *
822  (1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff)));
823  }
824  break;
825 
826  default:
827  break;
828  }
829  }
830 
831  return outarray;
832 }
833 
834 void MMFDiffusion::v_EvaluateExactSolution(unsigned int field,
835  Array<OneD, NekDouble> &outfield,
836  const NekDouble time)
837 {
838  switch (m_TestType)
839  {
840  case eTestPlane:
841  {
842  TestPlaneProblem(time, outfield);
843  }
844  break;
845 
846  case eTestCube:
847  {
848  TestCubeProblem(time, outfield);
849  }
850  break;
851 
852  case eTestLinearSphere:
853  case eTestNonlinearSphere:
854  {
855  Morphogenesis(time, field, outfield);
856  }
857  break;
858 
859  case eFHNStandard:
860  case eFHNRogers:
861  case eFHNAlievPanf:
862  {
863  int nq = GetTotPoints();
864  outfield = Array<OneD, NekDouble>(nq, 0.0);
865  }
866  /* Falls through. */
867  default:
868  {
869  EquationSystem::v_EvaluateExactSolution(field, outfield, time);
870  }
871  break;
872  }
873 }
874 
875 void MMFDiffusion::v_GenerateSummary(SolverUtils::SummaryList &s)
876 {
877  MMFSystem::v_GenerateSummary(s);
878  SolverUtils::AddSummaryItem(s, "TestType", TestTypeMap[m_TestType]);
879  SolverUtils::AddSummaryItem(s, "epsilon0", m_epsilon[0]);
880  SolverUtils::AddSummaryItem(s, "epsilon1", m_epsilon[1]);
881  SolverUtils::AddSummaryItem(s, "epsilon2", m_epsilon[2]);
882  if (m_TestType == eTestLinearSphere)
883  {
884  SolverUtils::AddSummaryItem(s, "epsilon for u", m_epsu[0]);
885  SolverUtils::AddSummaryItem(s, "epsilon for v", m_epsu[1]);
886  }
887 }
888 } // namespace Nektar
889 int main(int argc, char *argv[])
890 {
891  LibUtilities::SessionReaderSharedPtr session;
892  SpatialDomains::MeshGraphSharedPtr graph;
893  std::string vDriverModule;
894  DriverSharedPtr drv;
895 
896  try
897  {
898  // Create session reader.
899  session = LibUtilities::SessionReader::CreateInstance(argc, argv);
900 
901  // Create MeshGraph
902  graph = SpatialDomains::MeshGraph::Read(session);
903 
904  // Create driver
905  session->LoadSolverInfo("Driver", vDriverModule, "Standard");
906  drv = GetDriverFactory().CreateInstance(vDriverModule, session, graph);
907 
908  // Execute driver
909  drv->Execute();
910 
911  // Finalise session
912  session->Finalise();
913  }
914  catch (const std::runtime_error &e)
915  {
916  return 1;
917  }
918  catch (const std::string &eStr)
919  {
920  std::cout << "Error: " << eStr << std::endl;
921  }
922 
923  return 0;
924 }
m_mu
NekDouble m_mu
Definition: CompressibleBL.cpp:86
main
int main(int argc, char *argv[])
Definition: MMFDiffusion.cpp:889
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:88
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineImplicitSolve
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:120
Nektar::MMFDiffusion::m_varcoeff
StdRegions::VarCoeffMap m_varcoeff
Variable diffusivity.
Definition: MMFDiffusion.h:150
Nektar::MMFDiffusion::~MMFDiffusion
virtual ~MMFDiffusion()
Desctructor.
Definition: MMFDiffusion.cpp:184
Nektar::MMFDiffusion::Morphogenesis
void Morphogenesis(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFDiffusion.cpp:569
Nektar::MMFDiffusion::m_InitWaveType
InitWaveType m_InitWaveType
Definition: MMFDiffusion.h:108
Nektar::MMFDiffusion::m_epsu
Array< OneD, NekDouble > m_epsu
Definition: MMFDiffusion.h:153
Nektar::MMFDiffusion::TestCubeProblem
void TestCubeProblem(const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFDiffusion.cpp:549
Nektar::MMFDiffusion::m_epsilon
Array< OneD, NekDouble > m_epsilon
Definition: MMFDiffusion.h:152
Nektar::MMFDiffusion::v_GenerateSummary
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Prints a summary of the model parameters.
Definition: MMFDiffusion.cpp:875
Nektar::MMFDiffusion::DoOdeRhs
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, const NekDouble time)
Computes the reaction terms and .
Definition: MMFDiffusion.cpp:249
Nektar::MMFDiffusion::PlanePhiWave
Array< OneD, NekDouble > PlanePhiWave()
Definition: MMFDiffusion.cpp:706
Nektar::MMFDiffusion::DoImplicitSolve
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray, NekDouble time, NekDouble lambda)
Solve for the diffusion term.
Definition: MMFDiffusion.cpp:194
Nektar::MMFDiffusion::TestPlaneProblem
void TestPlaneProblem(const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFDiffusion.cpp:529
Nektar::MMFDiffusion::v_EvaluateExactSolution
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Definition: MMFDiffusion.cpp:834
Nektar::MMFDiffusion::v_SetInitialConditions
virtual void v_SetInitialConditions(NekDouble initialtime, bool dumpInitialConditions, const int domain)
Sets a custom initial condition.
Definition: MMFDiffusion.cpp:448
Nektar::MMFDiffusion::v_InitObject
virtual void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
Definition: MMFDiffusion.cpp:63
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:368
Nektar::SolverUtils::EquationSystem::v_SetInitialConditions
virtual SOLVER_UTILS_EXPORT void v_SetInitialConditions(NekDouble initialtime=0.0, bool dumpInitialConditions=true, const int domain=0)
Definition: EquationSystem.cpp:926
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:335
Nektar::SolverUtils::EquationSystem::v_EvaluateExactSolution
virtual SOLVER_UTILS_EXPORT void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time)
Definition: EquationSystem.cpp:991
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1123
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:341
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:328
Nektar::SolverUtils::EquationSystem::SetBoundaryConditions
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
Definition: EquationSystem.cpp:715
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:720
Nektar::SolverUtils::MMFSystem
A base class for PDEs which include an advection component.
Definition: MMFSystem.h:147
Nektar::SolverUtils::MMFSystem::m_DivMF
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:195
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:186
Nektar::SolverUtils::MMFSystem::MMFInitObject
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble >> &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:53
Nektar::SolverUtils::MMFSystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2469
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition: UnsteadySystem.h:49
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:85
Nektar::SolverUtils::UnsteadySystem::m_explicitDiffusion
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
Definition: UnsteadySystem.h:89
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true)
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:79
Nektar::LibUtilities::rad
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter.cpp:86
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::SolverUtils::DriverSharedPtr
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:51
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:48
Nektar::SolverUtils::GetDriverFactory
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:85
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::InitWaveTypeMap
const char *const InitWaveTypeMap[]
Definition: MMFDiffusion.h:74
Nektar::eLeftBottomCorner
@ eLeftBottomCorner
Definition: MMFDiffusion.h:68
Nektar::eBothEnds
@ eBothEnds
Definition: MMFDiffusion.h:66
Nektar::eSpiralDock
@ eSpiralDock
Definition: MMFDiffusion.h:70
Nektar::eFHNStandard
@ eFHNStandard
Definition: MMFDiffusion.h:52
Nektar::eTestLinearSphere
@ eTestLinearSphere
Definition: MMFDiffusion.h:50
Nektar::eTestPlane
@ eTestPlane
Definition: MMFDiffusion.h:48
Nektar::eTestCube
@ eTestCube
Definition: MMFDiffusion.h:49
Nektar::SIZE_TestType
@ SIZE_TestType
Length of enum list.
Definition: MMFDiffusion.h:55
Nektar::eTestNonlinearSphere
@ eTestNonlinearSphere
Definition: MMFDiffusion.h:51
Nektar::eFHNRogers
@ eFHNRogers
Definition: MMFDiffusion.h:53
Nektar::eFHNAlievPanf
@ eFHNAlievPanf
Definition: MMFDiffusion.h:54
Nektar::TestTypeMap
const char *const TestTypeMap[]
Definition: MMFDiffusion.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmul
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:209
Vmath::Svtvp
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:622
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:1050
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:574
Vmath::Svtvm
void Svtvm(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:664
Vmath::Vadd
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:359
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
Vmath::Vdiv
void Vdiv(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:284
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:384
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:945
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291
Generated on Tue Aug 23 2022 07:27:13 for Nektar++ by doxygen 1.9.1