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