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