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