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  Vmath::Vcopy(nq, &m_movingframes[k][j * nq], 1, &tmp[0], 1);
152  m_varcoeff[MMFCoeffs[indx + j]] = tmp;
153  }
154 
155  // m_DivMF
156  Vmath::Vcopy(nq, &m_DivMF[k][0], 1, &tmp[0], 1);
157  m_varcoeff[MMFCoeffs[indx + 3]] = tmp;
158 
159  // \| e^k \|
160  tmp = Array<OneD, NekDouble>(nq, 0.0);
161  for (int i = 0; i < m_spacedim; ++i)
162  {
163  Vmath::Vvtvp(nq, &m_movingframes[k][i * nq], 1,
164  &m_movingframes[k][i * nq], 1, &tmp[0], 1, &tmp[0], 1);
165  }
166 
167  m_varcoeff[MMFCoeffs[indx + 4]] = tmp;
168  }
169 
170  if (!m_explicitDiffusion)
171  {
172  m_ode.DefineImplicitSolve(&MMFDiffusion::DoImplicitSolve, this);
173  }
174 
175  m_ode.DefineOdeRhs(&MMFDiffusion::DoOdeRhs, this);
176 }
177 
178 /**
179  *
180  */
181 MMFDiffusion::~MMFDiffusion()
182 {
183 }
184 
185 /**OdeRhs
186  * @param inarray Input array.
187  * @param outarray Output array.
188  * @param time Current simulation time.
189  * @param lambda Timestep.
190  */
191 void MMFDiffusion::DoImplicitSolve(
192  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
193  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time,
194  const NekDouble lambda)
195 {
196  int nvariables = inarray.size();
197  int nq = m_fields[0]->GetNpoints();
198 
199  StdRegions::ConstFactorMap factors;
200  factors[StdRegions::eFactorTau] = 1.0;
201 
202  Array<OneD, Array<OneD, NekDouble>> F(nvariables);
203  factors[StdRegions::eFactorLambda] = 1.0 / lambda;
204  F[0] = Array<OneD, NekDouble>(nq * nvariables);
205 
206  for (int n = 1; n < nvariables; ++n)
207  {
208  F[n] = F[n - 1] + nq;
209  // cout << "F["<< n<<"=" << F[n][1] <<endl;
210  }
211 
212  // We solve ( \nabla^2 - HHlambda ) Y[i] = rhs [i]
213  // inarray = input: \hat{rhs} -> output: \hat{Y}
214  // outarray = output: nabla^2 \hat{Y}
215  // where \hat = modal coeffs
216  SetBoundaryConditions(time);
217 
218  for (int i = 0; i < nvariables; ++i)
219  {
220  factors[StdRegions::eFactorLambda] = 1.0 / lambda / m_epsu[i];
221 
222  // Multiply 1.0/timestep
223  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda], inarray[i], 1,
224  F[i], 1);
225 
226  /* for (int k = 0; k < 15; ++k)
227  cout << "inarray["<<i << "]"<< k<<"=" << inarray[i][k]<<endl;*/
228  // Solve a system of equations with Helmholtz solver and transform
229  // back into physical space.
230  m_fields[i]->HelmSolve(F[i], m_fields[i]->UpdateCoeffs(), factors,
231  m_varcoeff);
232 
233  m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
234  /* Array<OneD, NekDouble> coefarray = m_fields[i]->GetCoeffs();
235  for (int k = 0; k < 15; ++k)
236  cout << "inarray["<< k<<"=" << coefarray[k]<<endl;*/
237  }
238  /* for (int kk = 0; kk < 15; ++kk)
239  cout << "inarray["<< kk<<"=" <<
240  m_varcoeff[StdRegions::eVarCoeffMF3Mag][kk]<<endl;*/
241 }
242 
243 /**
244  *
245  */
246 void MMFDiffusion::DoOdeRhs(
247  const Array<OneD, const Array<OneD, NekDouble>> &inarray,
248  Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
249 {
250  int nq = GetTotPoints();
251 
252  switch (m_TestType)
253  {
254  case eTestPlane:
255  {
256 
257  Array<OneD, NekDouble> x(nq);
258  Array<OneD, NekDouble> y(nq);
259  Array<OneD, NekDouble> z(nq);
260 
261  m_fields[0]->GetCoords(x, y, z);
262 
263  for (int k = 0; k < nq; k++)
264  {
265  outarray[0][k] = (m_epsilon[0] + m_epsilon[1] - 1.0) * m_pi *
266  m_pi * exp(-1.0 * m_pi * m_pi * time) *
267  sin(m_pi * x[k]) * cos(m_pi * y[k]);
268  }
269  }
270  break;
271 
272  case eTestCube:
273  {
274 
275  Array<OneD, NekDouble> x(nq);
276  Array<OneD, NekDouble> y(nq);
277  Array<OneD, NekDouble> z(nq);
278 
279  m_fields[0]->GetCoords(x, y, z);
280 
281  for (int k = 0; k < nq; k++)
282  {
283  outarray[0][k] =
284  (m_epsilon[0] + m_epsilon[1] + m_epsilon[2] - 1.0) * m_pi *
285  m_pi * exp(-1.0 * m_pi * m_pi * time) * sin(m_pi * x[k]) *
286  sin(m_pi * y[k]) * sin(m_pi * z[k]);
287  }
288  }
289  break;
290 
291  case eTestLinearSphere:
292  {
293  Array<OneD, NekDouble> temp(nq);
294 
295  NekDouble A = 2.0;
296  NekDouble B = 5.0;
297 
298  NekDouble m_a, m_b, m_c, m_d;
299  m_a = B - 1.0;
300  m_b = A * A;
301  m_c = -1.0 * B;
302  m_d = -1.0 * A * A;
303 
304  temp = Array<OneD, NekDouble>(nq, 0.0);
305  Vmath::Svtvp(nq, m_a, &inarray[0][0], 1, &temp[0], 1, &temp[0], 1);
306  Vmath::Svtvp(nq, m_b, &inarray[1][0], 1, &temp[0], 1,
307  &outarray[0][0], 1);
308 
309  temp = Array<OneD, NekDouble>(nq, 0.0);
310  Vmath::Svtvp(nq, m_c, &inarray[0][0], 1, &temp[0], 1, &temp[0], 1);
311  Vmath::Svtvp(nq, m_d, &inarray[1][0], 1, &temp[0], 1,
312  &outarray[1][0], 1);
313  }
314  break;
315 
316  case eTestNonlinearSphere:
317  {
318  NekDouble A = 2.0;
319  NekDouble B = 5.0;
320 
321  Array<OneD, NekDouble> Aonevec(nq, A);
322 
323  // cube = phys0*phys0*phy1
324  Array<OneD, NekDouble> cube(nq);
325  Vmath::Vmul(nq, &inarray[0][0], 1, &inarray[0][0], 1, &cube[0], 1);
326  Vmath::Vmul(nq, &inarray[1][0], 1, &cube[0], 1, &cube[0], 1);
327 
328  // outarray[0] = A - B*phy0 + phy0*phy0*phy1 - phy0
329  NekDouble coeff = -1.0 * B - 1.0;
330  Array<OneD, NekDouble> tmp(nq);
331  Vmath::Svtvp(nq, coeff, &inarray[0][0], 1, &cube[0], 1, &tmp[0], 1);
332  Vmath::Vadd(nq, &Aonevec[0], 1, &tmp[0], 1, &outarray[0][0], 1);
333 
334  // outarray[1] = B*phys0 - phy0*phy0*phy1
335  Vmath::Svtvm(nq, B, &inarray[0][0], 1, &cube[0], 1, &outarray[1][0],
336  1);
337  }
338  break;
339 
340  case eFHNStandard:
341  {
342  // \phi - \phi^3/3 - \psi
343  NekDouble a = 0.12;
344  NekDouble b = 0.011;
345  NekDouble c1 = 0.175;
346  NekDouble c2 = 0.03;
347  NekDouble d = 0.55;
348 
349  Array<OneD, NekDouble> tmp(nq);
350 
351  // Reaction for \phi = c1 \phi ( \phi - a)*(1 - \phi) - c2 v
352  Vmath::Smul(nq, -1.0 * c1, inarray[0], 1, outarray[0], 1);
353  Vmath::Sadd(nq, -1.0 * a, inarray[0], 1, tmp, 1);
354  Vmath::Vmul(nq, tmp, 1, inarray[0], 1, outarray[0], 1);
355  Vmath::Sadd(nq, -1.0, inarray[0], 1, tmp, 1);
356  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
357 
358  Vmath::Smul(nq, -1.0 * c2, inarray[1], 1, tmp, 1);
359  Vmath::Vadd(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
360 
361  // Reaction for \psi = b (\phi - d \psi )
362  Vmath::Svtvp(nq, -1.0 * d, inarray[1], 1, inarray[0], 1,
363  outarray[1], 1);
364  Vmath::Smul(nq, b, outarray[1], 1, outarray[1], 1);
365  }
366  break;
367 
368  case eFHNRogers:
369  {
370  NekDouble a = 0.13;
371  NekDouble b = 0.013;
372  NekDouble c1 = 0.26;
373  NekDouble c2 = 0.1;
374  NekDouble d = 1.0;
375 
376  Array<OneD, NekDouble> tmp(nq);
377 
378  // Reaction for \phi = c1 \phi ( \phi - a)*(1 - \phi) - c2 u v
379  Vmath::Smul(nq, -1.0 * c1, inarray[0], 1, outarray[0], 1);
380  Vmath::Sadd(nq, -1.0 * a, inarray[0], 1, tmp, 1);
381  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
382  Vmath::Sadd(nq, -1.0, inarray[0], 1, tmp, 1);
383  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
384 
385  Vmath::Vmul(nq, inarray[0], 1, inarray[1], 1, tmp, 1);
386  Vmath::Smul(nq, -1.0 * c2, tmp, 1, tmp, 1);
387  Vmath::Vadd(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
388 
389  // Reaction for \psi = b (\phi - d \psi )
390  Vmath::Svtvp(nq, -1.0 * d, inarray[1], 1, inarray[0], 1,
391  outarray[1], 1);
392  Vmath::Smul(nq, b, outarray[1], 1, outarray[1], 1);
393  }
394  break;
395 
396  case eFHNAlievPanf:
397  {
398 
399  NekDouble a = 0.15;
400  NekDouble c1 = 8.0;
401  NekDouble c2 = 1.0;
402  NekDouble c0 = 0.002;
403  NekDouble mu1 = 0.2;
404  NekDouble mu2 = 0.3;
405 
406  Array<OneD, NekDouble> tmp(nq);
407 
408  // Reaction for \phi = c1 \phi ( \phi - a)*(1 - \phi) - c2 u v
409  Vmath::Smul(nq, -1.0 * c1, inarray[0], 1, outarray[0], 1);
410  Vmath::Sadd(nq, -1.0 * a, inarray[0], 1, tmp, 1);
411  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
412  Vmath::Sadd(nq, -1.0, inarray[0], 1, tmp, 1);
413  Vmath::Vmul(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
414 
415  Vmath::Vmul(nq, inarray[0], 1, inarray[1], 1, tmp, 1);
416  Vmath::Smul(nq, -1.0 * c2, tmp, 1, tmp, 1);
417  Vmath::Vadd(nq, tmp, 1, outarray[0], 1, outarray[0], 1);
418 
419  // Reaction for \psi = (c0 + (\mu1 \psi/(\mu2+\phi) ) )*(-\psi - c1
420  // * \phi*(\phi - a - 1) )
421 
422  Vmath::Smul(nq, mu1, inarray[1], 1, outarray[1], 1);
423  Vmath::Sadd(nq, mu2, inarray[0], 1, tmp, 1);
424  Vmath::Vdiv(nq, outarray[1], 1, tmp, 1, outarray[1], 1);
425  Vmath::Sadd(nq, c0, outarray[1], 1, outarray[1], 1);
426 
427  Vmath::Sadd(nq, (-a - 1.0), inarray[0], 1, tmp, 1);
428  Vmath::Vmul(nq, inarray[0], 1, tmp, 1, tmp, 1);
429  Vmath::Smul(nq, c1, tmp, 1, tmp, 1);
430  Vmath::Vadd(nq, inarray[1], 1, tmp, 1, tmp, 1);
431  Vmath::Neg(nq, tmp, 1);
432 
433  Vmath::Vmul(nq, tmp, 1, outarray[1], 1, outarray[1], 1);
434  }
435  break;
436 
437  default:
438  break;
439  }
440 }
441 
442 /**
443  *
444  */
445 void MMFDiffusion::v_SetInitialConditions(NekDouble initialtime,
446  bool dumpInitialConditions,
447  const int domain)
448 {
449  boost::ignore_unused(domain);
450 
451  int nq = GetTotPoints();
452 
453  switch (m_TestType)
454  {
455  case eTestPlane:
456  {
457  Array<OneD, NekDouble> u(nq);
458 
459  TestPlaneProblem(initialtime, u);
460  m_fields[0]->SetPhys(u);
461  }
462  break;
463 
464  case eTestCube:
465  {
466  Array<OneD, NekDouble> u(nq);
467 
468  TestCubeProblem(initialtime, u);
469  m_fields[0]->SetPhys(u);
470  /*for (int k=0; k<nq; ++k)
471  {
472  //for (int j=0; j<m_spacedim; ++j)
473  //{
474  cout << "_varcoeff" << u[k] <<endl;
475  // }
476  }*/
477  }
478  break;
479 
480  case eTestLinearSphere:
481  case eTestNonlinearSphere:
482  {
483  Array<OneD, NekDouble> u(nq);
484  Array<OneD, NekDouble> v(nq);
485 
486  Morphogenesis(initialtime, 0, u);
487  Morphogenesis(initialtime, 1, v);
488 
489  m_fields[0]->SetPhys(u);
490  m_fields[1]->SetPhys(v);
491  }
492  break;
493 
494  case eFHNStandard:
495  case eFHNRogers:
496  case eFHNAlievPanf:
497  {
498  Array<OneD, NekDouble> Zero(nq, 0.0);
499  m_fields[0]->SetPhys(PlanePhiWave());
500  m_fields[1]->SetPhys(Zero);
501  }
502  break;
503 
504  default:
505  {
506  EquationSystem::v_SetInitialConditions(initialtime, false);
507  }
508  break;
509  }
510 
511  // forward transform to fill the modal coeffs
512  for (int i = 0; i < m_fields.size(); ++i)
513  {
514  m_fields[i]->SetPhysState(true);
515  m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
516  m_fields[i]->UpdateCoeffs());
517  }
518 
519  if (dumpInitialConditions)
520  {
521  std::string outname = m_sessionName + "_initial.chk";
522  WriteFld(outname);
523  }
524 }
525 
526 void MMFDiffusion::TestPlaneProblem(const NekDouble time,
527  Array<OneD, NekDouble> &outfield)
528 
529 {
530  int nq = GetTotPoints();
531 
532  Array<OneD, NekDouble> x(nq);
533  Array<OneD, NekDouble> y(nq);
534  Array<OneD, NekDouble> z(nq);
535 
536  m_fields[0]->GetCoords(x, y, z);
537 
538  outfield = Array<OneD, NekDouble>(nq);
539  for (int k = 0; k < nq; k++)
540  {
541  outfield[k] = exp(-1.0 * m_pi * m_pi * time) * sin(m_pi * x[k]) *
542  cos(m_pi * y[k]);
543  }
544 }
545 
546 void MMFDiffusion::TestCubeProblem(const NekDouble time,
547  Array<OneD, NekDouble> &outfield)
548 
549 {
550  int nq = GetTotPoints();
551 
552  Array<OneD, NekDouble> x(nq);
553  Array<OneD, NekDouble> y(nq);
554  Array<OneD, NekDouble> z(nq);
555 
556  m_fields[0]->GetCoords(x, y, z);
557 
558  outfield = Array<OneD, NekDouble>(nq);
559  for (int k = 0; k < nq; k++)
560  {
561  outfield[k] = exp(-1.0 * m_pi * m_pi * time) * sin(m_pi * x[k]) *
562  sin(m_pi * y[k]) * sin(m_pi * z[k]);
563  }
564 }
565 
566 void MMFDiffusion::Morphogenesis(const NekDouble time, unsigned int field,
567  Array<OneD, NekDouble> &outfield)
568 {
569  int nq = GetTotPoints();
570 
571  int i, m, n, ind;
572  NekDouble a_n, d_n, gamma_n;
573  NekDouble A_mn, C_mn, theta, phi, radius;
574 
575  std::complex<double> Spericharmonic, delta_n, temp;
576  std::complex<double> varphi0, varphi1;
577  std::complex<double> B_mn, D_mn;
578 
579  // Set some parameter values
580  int Maxn = 6;
581  int Maxm = 2 * Maxn - 1;
582 
583  NekDouble A = 2.0;
584  NekDouble B = 5.0;
585 
586  NekDouble m_mu = 0.001;
587  NekDouble m_nu = 0.002;
588 
589  NekDouble m_a, m_b, m_c, m_d;
590 
591  m_a = B - 1.0;
592  m_b = A * A;
593  m_c = -1.0 * B;
594  m_d = -1.0 * A * A;
595 
596  Array<OneD, Array<OneD, NekDouble>> Ainit(Maxn);
597  Array<OneD, Array<OneD, NekDouble>> Binit(Maxn);
598 
599  for (i = 0; i < Maxn; ++i)
600  {
601  Ainit[i] = Array<OneD, NekDouble>(Maxm, 0.0);
602  Binit[i] = Array<OneD, NekDouble>(Maxm, 0.0);
603  }
604 
605  Ainit[5][0] = -0.5839;
606  Ainit[5][1] = -0.8436;
607  Ainit[5][2] = -0.4764;
608  Ainit[5][3] = 0.6475;
609  Ainit[5][4] = 0.1886;
610  Ainit[5][5] = 0.8709;
611  Ainit[5][6] = -0.8338;
612  Ainit[5][7] = 0.1795;
613  Ainit[5][8] = -0.7873;
614  Ainit[5][9] = 0.8842;
615  Ainit[5][10] = 0.2943;
616 
617  Binit[5][0] = -0.6263;
618  Binit[5][1] = 0.9803;
619  Binit[5][2] = 0.7222;
620  Binit[5][3] = 0.5945;
621  Binit[5][4] = 0.6026;
622  Binit[5][5] = -0.2076;
623  Binit[5][6] = 0.4556;
624  Binit[5][7] = 0.6024;
625  Binit[5][8] = 0.9695;
626  Binit[5][9] = -0.4936;
627  Binit[5][10] = 0.1098;
628 
629  Array<OneD, NekDouble> u(nq);
630  Array<OneD, NekDouble> v(nq);
631  Array<OneD, NekDouble> x(nq);
632  Array<OneD, NekDouble> y(nq);
633  Array<OneD, NekDouble> z(nq);
634 
635  m_fields[0]->GetCoords(x, y, z);
636  for (int i = 0; i < nq; ++i)
637  {
638  radius = sqrt(x[i] * x[i] + y[i] * y[i] + z[i] * z[i]);
639 
640  // theta is in [0, pi]
641  theta = asin(z[i] / radius) + 0.5 * m_pi;
642 
643  // phi is in [0, 2*pi]
644  phi = atan2(y[i], x[i]) + m_pi;
645 
646  varphi0 = 0.0 * varphi0;
647  varphi1 = 0.0 * varphi1;
648  for (n = 0; n < Maxn; ++n)
649  {
650  // Set up parameters
651  a_n = m_a - m_mu * (n * (n + 1) / radius / radius);
652  d_n = m_d - m_nu * (n * (n + 1) / radius / radius);
653 
654  gamma_n = 0.5 * (a_n + d_n);
655 
656  temp = (a_n + d_n) * (a_n + d_n) - 4.0 * (a_n * d_n - m_b * m_c);
657  delta_n = 0.5 * sqrt(temp);
658 
659  for (m = -n; m <= n; ++m)
660  {
661  ind = m + n;
662  A_mn = Ainit[n][ind];
663  C_mn = Binit[n][ind];
664 
665  B_mn = ((a_n - gamma_n) * Ainit[n][ind] + m_b * Binit[n][ind]) /
666  delta_n;
667  D_mn = (m_c * Ainit[n][ind] + (d_n - gamma_n) * Binit[n][ind]) /
668  delta_n;
669 
670  Spericharmonic =
671  boost::math::spherical_harmonic(n, m, theta, phi);
672  varphi0 += exp(gamma_n * time) *
673  (A_mn * cosh(delta_n * time) +
674  B_mn * sinh(delta_n * time)) *
675  Spericharmonic;
676  varphi1 += exp(gamma_n * time) *
677  (C_mn * cosh(delta_n * time) +
678  D_mn * sinh(delta_n * time)) *
679  Spericharmonic;
680  }
681  }
682 
683  u[i] = varphi0.real();
684  v[i] = varphi1.real();
685  }
686 
687  switch (field)
688  {
689  case 0:
690  {
691  outfield = u;
692  }
693  break;
694 
695  case 1:
696  {
697  outfield = v;
698  }
699  break;
700  }
701 }
702 
703 Array<OneD, NekDouble> MMFDiffusion::PlanePhiWave()
704 {
705  int nq = GetTotPoints();
706  Array<OneD, NekDouble> outarray(nq, 0.0);
707 
708  Array<OneD, NekDouble> x(nq);
709  Array<OneD, NekDouble> y(nq);
710  Array<OneD, NekDouble> z(nq);
711 
712  m_fields[0]->GetCoords(x, y, z);
713 
714  NekDouble xmin, ymin, xmax;
715 
716  xmin = Vmath::Vmin(nq, x, 1);
717  xmax = Vmath::Vmax(nq, x, 1);
718  ymin = Vmath::Vmin(nq, y, 1);
719 
720  NekDouble xp, yp, xp2;
721  for (int i = 0; i < nq; i++)
722  {
723  switch (m_InitWaveType)
724  {
725  case eLeft:
726  {
727  NekDouble radiusofinit = 4.0;
728  NekDouble frontstiff = 0.1;
729 
730  xp = x[i] - xmin;
731  outarray[i] =
732  1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff));
733  }
734  break;
735 
736  case eBothEnds:
737  {
738  NekDouble radiusofinit = 3.0;
739  NekDouble frontstiff = 0.1;
740 
741  xp = x[i] - xmin;
742  xp2 = x[i] - xmax;
743 
744  outarray[i] =
745  1.0 / (1.0 +
746  exp((sqrt(xp * xp) - radiusofinit) / frontstiff)) +
747  1.0 / (1.0 +
748  exp((sqrt(xp2 * xp2) - radiusofinit) / frontstiff));
749  }
750  break;
751 
752  case eCenter:
753  {
754  NekDouble radiusofinit = 6.0;
755  NekDouble frontstiff = 0.1;
756 
757  // NekDouble xc = 0.5*(Vmath::Vmax(nq, x, 1) + Vmath::Vmin(nq,
758  // x, 1));
759 
760  xp = x[i] - xmin;
761  outarray[i] =
762  1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff));
763  }
764  break;
765 
766  case eLeftBottomCorner:
767  {
768  NekDouble radiusofinit = 6.0;
769  NekDouble frontstiff = 0.1;
770  NekDouble bs = 2.0;
771 
772  xp = x[i] - xmin;
773  yp = y[i] - ymin;
774  outarray[i] =
775  1.0 /
776  (1.0 + exp((sqrt(xp * xp + yp * yp) / bs - radiusofinit) /
777  frontstiff));
778  }
779  break;
780 
781  case ePoint:
782  {
783  NekDouble xloc, yloc, zloc, rad;
784  NekDouble radiusofinit = 10.0;
785 
786  xloc = x[i] - m_InitPtx;
787  yloc = y[i] - m_InitPty;
788  zloc = z[i] - m_InitPtz;
789 
790  rad = sqrt(xloc * xloc + yloc * yloc + zloc * zloc);
791 
792  xloc = xloc / radiusofinit;
793  yloc = yloc / radiusofinit;
794  zloc = zloc / radiusofinit;
795 
796  if (rad < radiusofinit)
797  {
798  outarray[i] =
799  exp(-(1.0 / 2.0) *
800  (xloc * xloc + yloc * yloc + zloc * zloc));
801  }
802 
803  else
804  {
805  outarray[i] = 0.0;
806  }
807  }
808  break;
809 
810  case eSpiralDock:
811  {
812  NekDouble radiusofinit = 3.0;
813  NekDouble frontstiff = 0.1;
814  xp = x[i] - 4.0;
815  yp = y[i];
816  outarray[i] =
817  (1.0 / (1.0 + exp(2.0 * yp))) *
818  (1.0 / (1.0 + exp(-2.0 * xp))) *
819  (1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff)));
820  }
821  break;
822 
823  default:
824  break;
825  }
826  }
827 
828  return outarray;
829 }
830 
831 void MMFDiffusion::v_EvaluateExactSolution(unsigned int field,
832  Array<OneD, NekDouble> &outfield,
833  const NekDouble time)
834 {
835  switch (m_TestType)
836  {
837  case eTestPlane:
838  {
839  TestPlaneProblem(time, outfield);
840  }
841  break;
842 
843  case eTestCube:
844  {
845  TestCubeProblem(time, outfield);
846  }
847  break;
848 
849  case eTestLinearSphere:
850  case eTestNonlinearSphere:
851  {
852  Morphogenesis(time, field, outfield);
853  }
854  break;
855 
856  case eFHNStandard:
857  case eFHNRogers:
858  case eFHNAlievPanf:
859  {
860  int nq = GetTotPoints();
861  outfield = Array<OneD, NekDouble>(nq, 0.0);
862  }
863  /* Falls through. */
864  default:
865  {
866  EquationSystem::v_EvaluateExactSolution(field, outfield, time);
867  }
868  break;
869  }
870 }
871 
872 void MMFDiffusion::v_GenerateSummary(SolverUtils::SummaryList &s)
873 {
874  MMFSystem::v_GenerateSummary(s);
875  SolverUtils::AddSummaryItem(s, "TestType", TestTypeMap[m_TestType]);
876  SolverUtils::AddSummaryItem(s, "epsilon0", m_epsilon[0]);
877  SolverUtils::AddSummaryItem(s, "epsilon1", m_epsilon[1]);
878  SolverUtils::AddSummaryItem(s, "epsilon2", m_epsilon[2]);
879  if (m_TestType == eTestLinearSphere)
880  {
881  SolverUtils::AddSummaryItem(s, "epsilon for u", m_epsu[0]);
882  SolverUtils::AddSummaryItem(s, "epsilon for v", m_epsu[1]);
883  }
884 }
885 } // namespace Nektar
886 int main(int argc, char *argv[])
887 {
888  LibUtilities::SessionReaderSharedPtr session;
889  SpatialDomains::MeshGraphSharedPtr graph;
890  std::string vDriverModule;
891  DriverSharedPtr drv;
892 
893  try
894  {
895  // Create session reader.
896  session = LibUtilities::SessionReader::CreateInstance(argc, argv);
897 
898  // Create MeshGraph
899  graph = SpatialDomains::MeshGraph::Read(session);
900 
901  // Create driver
902  session->LoadSolverInfo("Driver", vDriverModule, "Standard");
903  drv = GetDriverFactory().CreateInstance(vDriverModule, session, graph);
904 
905  // Execute driver
906  drv->Execute();
907 
908  // Finalise session
909  session->Finalise();
910  }
911  catch (const std::runtime_error &e)
912  {
913  return 1;
914  }
915  catch (const std::string &eStr)
916  {
917  std::cout << "Error: " << eStr << std::endl;
918  }
919 
920  return 0;
921 }
m_mu
NekDouble m_mu
Definition: CompressibleBL.cpp:86
main
int main(int argc, char *argv[])
Definition: MMFDiffusion.cpp:886
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:104
Nektar::MMFDiffusion::m_varcoeff
StdRegions::VarCoeffMap m_varcoeff
Variable diffusivity.
Definition: MMFDiffusion.h:150
Nektar::MMFDiffusion::~MMFDiffusion
virtual ~MMFDiffusion()
Desctructor.
Definition: MMFDiffusion.cpp:181
Nektar::MMFDiffusion::v_GenerateSummary
virtual void v_GenerateSummary(SolverUtils::SummaryList &s) override
Prints a summary of the model parameters.
Definition: MMFDiffusion.cpp:872
Nektar::MMFDiffusion::Morphogenesis
void Morphogenesis(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFDiffusion.cpp:566
Nektar::MMFDiffusion::v_InitObject
virtual void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: MMFDiffusion.cpp:63
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:546
Nektar::MMFDiffusion::m_epsilon
Array< OneD, NekDouble > m_epsilon
Definition: MMFDiffusion.h:152
Nektar::MMFDiffusion::v_SetInitialConditions
virtual void v_SetInitialConditions(NekDouble initialtime, bool dumpInitialConditions, const int domain) override
Sets a custom initial condition.
Definition: MMFDiffusion.cpp:445
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:246
Nektar::MMFDiffusion::PlanePhiWave
Array< OneD, NekDouble > PlanePhiWave()
Definition: MMFDiffusion.cpp:703
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:191
Nektar::MMFDiffusion::TestPlaneProblem
void TestPlaneProblem(const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFDiffusion.cpp:526
Nektar::MMFDiffusion::v_EvaluateExactSolution
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Definition: MMFDiffusion.cpp:831
Nektar::SolverUtils::EquationSystem::m_spacedim
int m_spacedim
Spatial dimension (>= expansion dim).
Definition: EquationSystem.h:412
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:962
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:375
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:1045
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1230
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:381
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:368
Nektar::SolverUtils::EquationSystem::SetBoundaryConditions
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
Definition: EquationSystem.cpp:751
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:768
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::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2469
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::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:83
Nektar::SolverUtils::UnsteadySystem::m_explicitDiffusion
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
Definition: UnsteadySystem.h:87
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:79
Nektar::LibUtilities::rad
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter/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:399
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
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 minus 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 scalar 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:294
Generated on Sat Mar 25 2023 12:35:51 for Nektar++ by doxygen 1.9.1