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