Nektar++
MMFDiffusion.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: MMFDiffusion.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: MMFDiffusion.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36#include <iostream>
37
38#include <boost/algorithm/string.hpp>
39#include <boost/core/ignore_unused.hpp>
40
41#include <DiffusionSolver/EquationSystems/MMFDiffusion.h>
42#include <LibUtilities/BasicUtils/SessionReader.h>
43#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
44#include <SolverUtils/Driver.h>
45
46#include <boost/math/special_functions/spherical_harmonic.hpp>
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 const int domain)
439{
440 boost::ignore_unused(domain);
441
442 int nq = GetTotPoints();
443
444 switch (m_TestType)
445 {
446 case eTestPlane:
447 {
448 Array<OneD, NekDouble> u(nq);
449
450 TestPlaneProblem(initialtime, u);
451 m_fields[0]->SetPhys(u);
452 }
453 break;
454
455 case eTestCube:
456 {
457 Array<OneD, NekDouble> u(nq);
458
459 TestCubeProblem(initialtime, u);
460 m_fields[0]->SetPhys(u);
461 }
462 break;
463
464 case eTestLinearSphere:
465 case eTestNonlinearSphere:
466 {
467 Array<OneD, NekDouble> u(nq);
468 Array<OneD, NekDouble> v(nq);
469
470 Morphogenesis(initialtime, 0, u);
471 Morphogenesis(initialtime, 1, v);
472
473 m_fields[0]->SetPhys(u);
474 m_fields[1]->SetPhys(v);
475 }
476 break;
477
478 case eFHNStandard:
479 case eFHNRogers:
480 case eFHNAlievPanf:
481 {
482 Array<OneD, NekDouble> Zero(nq, 0.0);
483 m_fields[0]->SetPhys(PlanePhiWave());
484 m_fields[1]->SetPhys(Zero);
485 }
486 break;
487
488 default:
489 {
490 EquationSystem::v_SetInitialConditions(initialtime, false);
491 }
492 break;
493 }
494
495 // forward transform to fill the modal coeffs
496 for (int i = 0; i < m_fields.size(); ++i)
497 {
498 m_fields[i]->SetPhysState(true);
499 m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
500 m_fields[i]->UpdateCoeffs());
501 }
502
503 if (dumpInitialConditions)
504 {
505 std::string outname = m_sessionName + "_initial.chk";
506 WriteFld(outname);
507 }
508}
509
510void MMFDiffusion::TestPlaneProblem(const NekDouble time,
511 Array<OneD, NekDouble> &outfield)
512
513{
514 int nq = GetTotPoints();
515
516 Array<OneD, NekDouble> x(nq);
517 Array<OneD, NekDouble> y(nq);
518 Array<OneD, NekDouble> z(nq);
519
520 m_fields[0]->GetCoords(x, y, z);
521
522 outfield = Array<OneD, NekDouble>(nq);
523 for (int k = 0; k < nq; k++)
524 {
525 outfield[k] = exp(-1.0 * m_pi * m_pi * time) * sin(m_pi * x[k]) *
526 cos(m_pi * y[k]);
527 }
528}
529
530void MMFDiffusion::TestCubeProblem(const NekDouble time,
531 Array<OneD, NekDouble> &outfield)
532
533{
534 int nq = GetTotPoints();
535
536 Array<OneD, NekDouble> x(nq);
537 Array<OneD, NekDouble> y(nq);
538 Array<OneD, NekDouble> z(nq);
539
540 m_fields[0]->GetCoords(x, y, z);
541
542 outfield = Array<OneD, NekDouble>(nq);
543 for (int k = 0; k < nq; k++)
544 {
545 outfield[k] = exp(-1.0 * m_pi * m_pi * time) * sin(m_pi * x[k]) *
546 sin(m_pi * y[k]) * sin(m_pi * z[k]);
547 }
548}
549
550void MMFDiffusion::Morphogenesis(const NekDouble time, unsigned int field,
551 Array<OneD, NekDouble> &outfield)
552{
553 int nq = GetTotPoints();
554
555 int i, m, n, ind;
556 NekDouble a_n, d_n, gamma_n;
557 NekDouble A_mn, C_mn, theta, phi, radius;
558
559 std::complex<double> Spericharmonic, delta_n, temp;
560 std::complex<double> varphi0, varphi1;
561 std::complex<double> B_mn, D_mn;
562
563 // Set some parameter values
564 int Maxn = 6;
565 int Maxm = 2 * Maxn - 1;
566
567 NekDouble A = 2.0;
568 NekDouble B = 5.0;
569
570 NekDouble m_mu = 0.001;
571 NekDouble m_nu = 0.002;
572
573 NekDouble m_a, m_b, m_c, m_d;
574
575 m_a = B - 1.0;
576 m_b = A * A;
577 m_c = -1.0 * B;
578 m_d = -1.0 * A * A;
579
580 Array<OneD, Array<OneD, NekDouble>> Ainit(Maxn);
581 Array<OneD, Array<OneD, NekDouble>> Binit(Maxn);
582
583 for (i = 0; i < Maxn; ++i)
584 {
585 Ainit[i] = Array<OneD, NekDouble>(Maxm, 0.0);
586 Binit[i] = Array<OneD, NekDouble>(Maxm, 0.0);
587 }
588
589 Ainit[5][0] = -0.5839;
590 Ainit[5][1] = -0.8436;
591 Ainit[5][2] = -0.4764;
592 Ainit[5][3] = 0.6475;
593 Ainit[5][4] = 0.1886;
594 Ainit[5][5] = 0.8709;
595 Ainit[5][6] = -0.8338;
596 Ainit[5][7] = 0.1795;
597 Ainit[5][8] = -0.7873;
598 Ainit[5][9] = 0.8842;
599 Ainit[5][10] = 0.2943;
600
601 Binit[5][0] = -0.6263;
602 Binit[5][1] = 0.9803;
603 Binit[5][2] = 0.7222;
604 Binit[5][3] = 0.5945;
605 Binit[5][4] = 0.6026;
606 Binit[5][5] = -0.2076;
607 Binit[5][6] = 0.4556;
608 Binit[5][7] = 0.6024;
609 Binit[5][8] = 0.9695;
610 Binit[5][9] = -0.4936;
611 Binit[5][10] = 0.1098;
612
613 Array<OneD, NekDouble> u(nq);
614 Array<OneD, NekDouble> v(nq);
615 Array<OneD, NekDouble> x(nq);
616 Array<OneD, NekDouble> y(nq);
617 Array<OneD, NekDouble> z(nq);
618
619 m_fields[0]->GetCoords(x, y, z);
620 for (int i = 0; i < nq; ++i)
621 {
622 radius = sqrt(x[i] * x[i] + y[i] * y[i] + z[i] * z[i]);
623
624 // theta is in [0, pi]
625 theta = asin(z[i] / radius) + 0.5 * m_pi;
626
627 // phi is in [0, 2*pi]
628 phi = atan2(y[i], x[i]) + m_pi;
629
630 varphi0 = 0.0 * varphi0;
631 varphi1 = 0.0 * varphi1;
632 for (n = 0; n < Maxn; ++n)
633 {
634 // Set up parameters
635 a_n = m_a - m_mu * (n * (n + 1) / radius / radius);
636 d_n = m_d - m_nu * (n * (n + 1) / radius / radius);
637
638 gamma_n = 0.5 * (a_n + d_n);
639
640 temp = (a_n + d_n) * (a_n + d_n) - 4.0 * (a_n * d_n - m_b * m_c);
641 delta_n = 0.5 * sqrt(temp);
642
643 for (m = -n; m <= n; ++m)
644 {
645 ind = m + n;
646 A_mn = Ainit[n][ind];
647 C_mn = Binit[n][ind];
648
649 B_mn = ((a_n - gamma_n) * Ainit[n][ind] + m_b * Binit[n][ind]) /
650 delta_n;
651 D_mn = (m_c * Ainit[n][ind] + (d_n - gamma_n) * Binit[n][ind]) /
652 delta_n;
653
654 Spericharmonic =
655 boost::math::spherical_harmonic(n, m, theta, phi);
656 varphi0 += exp(gamma_n * time) *
657 (A_mn * cosh(delta_n * time) +
658 B_mn * sinh(delta_n * time)) *
659 Spericharmonic;
660 varphi1 += exp(gamma_n * time) *
661 (C_mn * cosh(delta_n * time) +
662 D_mn * sinh(delta_n * time)) *
663 Spericharmonic;
664 }
665 }
666
667 u[i] = varphi0.real();
668 v[i] = varphi1.real();
669 }
670
671 switch (field)
672 {
673 case 0:
674 {
675 outfield = u;
676 }
677 break;
678
679 case 1:
680 {
681 outfield = v;
682 }
683 break;
684 }
685}
686
687Array<OneD, NekDouble> MMFDiffusion::PlanePhiWave()
688{
689 int nq = GetTotPoints();
690 Array<OneD, NekDouble> outarray(nq, 0.0);
691
692 Array<OneD, NekDouble> x(nq);
693 Array<OneD, NekDouble> y(nq);
694 Array<OneD, NekDouble> z(nq);
695
696 m_fields[0]->GetCoords(x, y, z);
697
698 NekDouble xmin, ymin, xmax;
699
700 xmin = Vmath::Vmin(nq, x, 1);
701 xmax = Vmath::Vmax(nq, x, 1);
702 ymin = Vmath::Vmin(nq, y, 1);
703
704 NekDouble xp, yp, xp2;
705 for (int i = 0; i < nq; i++)
706 {
707 switch (m_InitWaveType)
708 {
709 case eLeft:
710 {
711 NekDouble radiusofinit = 4.0;
712 NekDouble frontstiff = 0.1;
713
714 xp = x[i] - xmin;
715 outarray[i] =
716 1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff));
717 }
718 break;
719
720 case eBothEnds:
721 {
722 NekDouble radiusofinit = 3.0;
723 NekDouble frontstiff = 0.1;
724
725 xp = x[i] - xmin;
726 xp2 = x[i] - xmax;
727
728 outarray[i] =
729 1.0 / (1.0 +
730 exp((sqrt(xp * xp) - radiusofinit) / frontstiff)) +
731 1.0 / (1.0 +
732 exp((sqrt(xp2 * xp2) - radiusofinit) / frontstiff));
733 }
734 break;
735
736 case eCenter:
737 {
738 NekDouble radiusofinit = 6.0;
739 NekDouble frontstiff = 0.1;
740
741 xp = x[i] - xmin;
742 outarray[i] =
743 1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff));
744 }
745 break;
746
747 case eLeftBottomCorner:
748 {
749 NekDouble radiusofinit = 6.0;
750 NekDouble frontstiff = 0.1;
751 NekDouble bs = 2.0;
752
753 xp = x[i] - xmin;
754 yp = y[i] - ymin;
755 outarray[i] =
756 1.0 /
757 (1.0 + exp((sqrt(xp * xp + yp * yp) / bs - radiusofinit) /
758 frontstiff));
759 }
760 break;
761
762 case ePoint:
763 {
764 NekDouble xloc, yloc, zloc, rad;
765 NekDouble radiusofinit = 10.0;
766
767 xloc = x[i] - m_InitPtx;
768 yloc = y[i] - m_InitPty;
769 zloc = z[i] - m_InitPtz;
770
771 rad = sqrt(xloc * xloc + yloc * yloc + zloc * zloc);
772
773 xloc = xloc / radiusofinit;
774 yloc = yloc / radiusofinit;
775 zloc = zloc / radiusofinit;
776
777 if (rad < radiusofinit)
778 {
779 outarray[i] =
780 exp(-(1.0 / 2.0) *
781 (xloc * xloc + yloc * yloc + zloc * zloc));
782 }
783
784 else
785 {
786 outarray[i] = 0.0;
787 }
788 }
789 break;
790
791 case eSpiralDock:
792 {
793 NekDouble radiusofinit = 3.0;
794 NekDouble frontstiff = 0.1;
795 xp = x[i] - 4.0;
796 yp = y[i];
797 outarray[i] =
798 (1.0 / (1.0 + exp(2.0 * yp))) *
799 (1.0 / (1.0 + exp(-2.0 * xp))) *
800 (1.0 / (1.0 + exp((xp - radiusofinit) / frontstiff)));
801 }
802 break;
803
804 default:
805 break;
806 }
807 }
808
809 return outarray;
810}
811
812void MMFDiffusion::v_EvaluateExactSolution(unsigned int field,
813 Array<OneD, NekDouble> &outfield,
814 const NekDouble time)
815{
816 switch (m_TestType)
817 {
818 case eTestPlane:
819 {
820 TestPlaneProblem(time, outfield);
821 }
822 break;
823
824 case eTestCube:
825 {
826 TestCubeProblem(time, outfield);
827 }
828 break;
829
830 case eTestLinearSphere:
831 case eTestNonlinearSphere:
832 {
833 Morphogenesis(time, field, outfield);
834 }
835 break;
836
837 case eFHNStandard:
838 case eFHNRogers:
839 case eFHNAlievPanf:
840 {
841 int nq = GetTotPoints();
842 outfield = Array<OneD, NekDouble>(nq, 0.0);
843 }
844 /* Falls through. */
845 default:
846 {
847 EquationSystem::v_EvaluateExactSolution(field, outfield, time);
848 }
849 break;
850 }
851}
852
853void MMFDiffusion::v_GenerateSummary(SolverUtils::SummaryList &s)
854{
855 MMFSystem::v_GenerateSummary(s);
856 SolverUtils::AddSummaryItem(s, "TestType", TestTypeMap[m_TestType]);
857 SolverUtils::AddSummaryItem(s, "epsilon0", m_epsilon[0]);
858 SolverUtils::AddSummaryItem(s, "epsilon1", m_epsilon[1]);
859 SolverUtils::AddSummaryItem(s, "epsilon2", m_epsilon[2]);
860 if (m_TestType == eTestLinearSphere)
861 {
862 SolverUtils::AddSummaryItem(s, "epsilon for u", m_epsu[0]);
863 SolverUtils::AddSummaryItem(s, "epsilon for v", m_epsu[1]);
864 }
865}
866} // namespace Nektar
867
868int main(int argc, char *argv[])
869{
870 LibUtilities::SessionReaderSharedPtr session;
871 SpatialDomains::MeshGraphSharedPtr graph;
872 std::string vDriverModule;
873 DriverSharedPtr drv;
874
875 try
876 {
877 // Create session reader.
878 session = LibUtilities::SessionReader::CreateInstance(argc, argv);
879
880 // Create MeshGraph
881 graph = SpatialDomains::MeshGraph::Read(session);
882
883 // Create driver
884 session->LoadSolverInfo("Driver", vDriverModule, "Standard");
885 drv = GetDriverFactory().CreateInstance(vDriverModule, session, graph);
886
887 // Execute driver
888 drv->Execute();
889
890 // Finalise session
891 session->Finalise();
892 }
893 catch (const std::runtime_error &e)
894 {
895 return 1;
896 }
897 catch (const std::string &eStr)
898 {
899 std::cout << "Error: " << eStr << std::endl;
900 }
901
902 return 0;
903}
m_mu
NekDouble m_mu
Definition: CompressibleBL.cpp:86
main
int main(int argc, char *argv[])
Definition: MMFDiffusion.cpp:868
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineOdeRhs
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:88
Nektar::LibUtilities::TimeIntegrationSchemeOperators::DefineImplicitSolve
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
Definition: TimeIntegrationSchemeOperators.h:104
Nektar::MMFDiffusion::m_varcoeff
StdRegions::VarCoeffMap m_varcoeff
Variable diffusivity.
Definition: MMFDiffusion.h:150
Nektar::MMFDiffusion::~MMFDiffusion
virtual ~MMFDiffusion()
Desctructor.
Definition: MMFDiffusion.cpp:181
Nektar::MMFDiffusion::v_GenerateSummary
virtual void v_GenerateSummary(SolverUtils::SummaryList &s) override
Prints a summary of the model parameters.
Definition: MMFDiffusion.cpp:853
Nektar::MMFDiffusion::Morphogenesis
void Morphogenesis(const NekDouble time, unsigned int field, Array< OneD, NekDouble > &outfield)
Definition: MMFDiffusion.cpp:550
Nektar::MMFDiffusion::v_InitObject
virtual void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: MMFDiffusion.cpp:63
Nektar::MMFDiffusion::m_InitWaveType
InitWaveType m_InitWaveType
Definition: MMFDiffusion.h:108
Nektar::MMFDiffusion::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:530
Nektar::MMFDiffusion::m_epsilon
Array< OneD, NekDouble > m_epsilon
Definition: MMFDiffusion.h:152
Nektar::MMFDiffusion::v_SetInitialConditions
virtual 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:687
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::TestPlaneProblem
void TestPlaneProblem(const NekDouble time, Array< OneD, NekDouble > &outfield)
Definition: MMFDiffusion.cpp:510
Nektar::MMFDiffusion::v_EvaluateExactSolution
virtual void v_EvaluateExactSolution(unsigned int field, Array< OneD, NekDouble > &outfield, const NekDouble time) override
Definition: MMFDiffusion.cpp:812
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:400
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:982
Nektar::SolverUtils::EquationSystem::m_fields
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
Definition: EquationSystem.h:363
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:1062
Nektar::SolverUtils::EquationSystem::WriteFld
SOLVER_UTILS_EXPORT void WriteFld(const std::string &outname)
Write field data to the given filename.
Definition: EquationSystem.cpp:1249
Nektar::SolverUtils::EquationSystem::m_sessionName
std::string m_sessionName
Name of the session.
Definition: EquationSystem.h:369
Nektar::SolverUtils::EquationSystem::m_session
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Definition: EquationSystem.h:356
Nektar::SolverUtils::EquationSystem::SetBoundaryConditions
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
Definition: EquationSystem.cpp:771
Nektar::SolverUtils::EquationSystem::GetTotPoints
SOLVER_UTILS_EXPORT int GetTotPoints()
Definition: EquationSystem.h:765
Nektar::SolverUtils::MMFSystem
A base class for PDEs which include an advection component.
Definition: MMFSystem.h:146
Nektar::SolverUtils::MMFSystem::m_DivMF
Array< OneD, Array< OneD, NekDouble > > m_DivMF
Definition: MMFSystem.h:194
Nektar::SolverUtils::MMFSystem::v_GenerateSummary
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
Definition: MMFSystem.cpp:2452
Nektar::SolverUtils::MMFSystem::m_movingframes
Array< OneD, Array< OneD, NekDouble > > m_movingframes
Definition: MMFSystem.h:185
Nektar::SolverUtils::MMFSystem::MMFInitObject
SOLVER_UTILS_EXPORT void MMFInitObject(const Array< OneD, const Array< OneD, NekDouble > > &Anisotropy, const int TangentXelem=-1)
Definition: MMFSystem.cpp:53
Nektar::SolverUtils::UnsteadySystem
Base class for unsteady solvers.
Definition: UnsteadySystem.h:48
Nektar::SolverUtils::UnsteadySystem::m_ode
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
Definition: UnsteadySystem.h:79
Nektar::SolverUtils::UnsteadySystem::m_explicitDiffusion
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
Definition: UnsteadySystem.h:97
Nektar::SolverUtils::UnsteadySystem::v_InitObject
virtual SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Init object for UnsteadySystem class.
Definition: UnsteadySystem.cpp:87
Nektar::LibUtilities::rad
static NekDouble rad(NekDouble x, NekDouble y)
Definition: Interpreter/Interpreter.cpp:86
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::SolverUtils::DriverSharedPtr
std::shared_ptr< Driver > DriverSharedPtr
A shared pointer to a Driver object.
Definition: Driver.h:54
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:48
Nektar::SolverUtils::GetDriverFactory
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:67
Nektar::SolverUtils::GetEquationSystemFactory
EquationSystemFactory & GetEquationSystemFactory()
Definition: EquationSystem.cpp:103
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:176
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:408
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
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
Nektar::InitWaveTypeMap
const char *const InitWaveTypeMap[]
Definition: MMFDiffusion.h:74
Nektar::eLeftBottomCorner
@ eLeftBottomCorner
Definition: MMFDiffusion.h:68
Nektar::eBothEnds
@ eBothEnds
Definition: MMFDiffusion.h:66
Nektar::eSpiralDock
@ eSpiralDock
Definition: MMFDiffusion.h:70
Nektar::eFHNStandard
@ eFHNStandard
Definition: MMFDiffusion.h:52
Nektar::eTestLinearSphere
@ eTestLinearSphere
Definition: MMFDiffusion.h:50
Nektar::eTestPlane
@ eTestPlane
Definition: MMFDiffusion.h:48
Nektar::eTestCube
@ eTestCube
Definition: MMFDiffusion.h:49
Nektar::SIZE_TestType
@ SIZE_TestType
Length of enum list.
Definition: MMFDiffusion.h:55
Nektar::eTestNonlinearSphere
@ eTestNonlinearSphere
Definition: MMFDiffusion.h:51
Nektar::eFHNRogers
@ eFHNRogers
Definition: MMFDiffusion.h:53
Nektar::eFHNAlievPanf
@ eFHNAlievPanf
Definition: MMFDiffusion.h:54
Nektar::TestTypeMap
const char *const TestTypeMap[]
Definition: MMFDiffusion.h:58
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:207
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:617
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:513
Vmath::Vmin
T Vmin(int n, const T *x, const int incx)
Return the minimum element in x - called vmin to avoid conflict with min.
Definition: Vmath.cpp:1045
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:569
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.cpp:659
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:354
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:245
Vmath::Vdiv
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
Definition: Vmath.cpp:280
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:487
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:43
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add scalar y = alpha + x.
Definition: Vmath.cpp:379
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:940
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1191
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Wed Jan 31 2024 08:26:50 for Nektar++ by doxygen 1.9.4