Skip to content

Commit c7a3191

Browse files
ProExpertProgProExpertProg
committed
Move kx/ky/kperp to grid, move hyper to own file
1 parent 0c8292c commit c7a3191

File tree

8 files changed

+120
-118
lines changed

8 files changed

+120
-118
lines changed

‎cpp/lib/Filter.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -6,16 +6,16 @@
66
namespace ahr {
77
void HouLiFilter::operator()(Grid::View::C_XY view) const {
88
grid.for_each_kxky([&](Dim kx, Dim ky) {
9-
view(kx, ky) *=
10-
exp(-36.0 * pow(kx_(kx) / grid.KX, 36.0)) * exp(-36.0 * pow(ky_(ky) / grid.KY, 36.0));
9+
view(kx, ky) *= exp(-36.0 * pow(grid.kx_(kx) / grid.KX, 36.0)) *
10+
exp(-36.0 * pow(grid.ky_(ky) / grid.KY, 36.0));
1111
});
1212
}
1313

1414
HouLiFilterCached::HouLiFilterCached(Grid const &grid)
1515
: HouLiFilter(grid), factors(std::array{grid.KX, grid.KY}) {
1616
grid.for_each_kxky([&](Dim kx, Dim ky) {
17-
factors(kx, ky) =
18-
exp(-36.0 * pow(kx_(kx) / grid.KX, 36.0)) * exp(-36.0 * pow(ky_(ky) / grid.KY, 36.0));
17+
factors(kx, ky) = exp(-36.0 * pow(grid.kx_(kx) / grid.KX, 36.0)) *
18+
exp(-36.0 * pow(grid.ky_(ky) / grid.KY, 36.0));
1919
});
2020
}
2121

@@ -26,10 +26,10 @@ void HouLiFilterCached::operator()(Grid::View::C_XY view) const {
2626
HouLiFilterCached1D::HouLiFilterCached1D(Grid const &grid)
2727
: HouLiFilter(grid), factors_x(grid.KX), factors_y(grid.KY) {
2828
for (Dim kx = 0; kx < grid.KX; ++kx) {
29-
factors_x.at(kx) = exp(-36.0 * pow(kx_(kx) / grid.KX, 36.0));
29+
factors_x.at(kx) = exp(-36.0 * pow(grid.kx_(kx) / grid.KX, 36.0));
3030
}
3131
for (Dim ky = 0; ky < grid.KY; ++ky) {
32-
factors_y.at(ky) = exp(-36.0 * pow(ky_(ky) / grid.KY, 36.0));
32+
factors_y.at(ky) = exp(-36.0 * pow(grid.ky_(ky) / grid.KY, 36.0));
3333
}
3434
}
3535

‎cpp/lib/Filter.hpp‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -12,12 +12,6 @@ class HouLiFilter {
1212

1313
protected:
1414
Grid const &grid;
15-
16-
// TODO extract to common utility
17-
[[nodiscard]] Real ky_(Dim ky) const {
18-
return (ky <= (grid.KY / 2) ? Real(ky) : Real(ky) - Real(grid.KY)) * Real(lx) / Real(ly);
19-
}
20-
[[nodiscard]] Real kx_(Dim kx) const { return Real(kx); }
2115
};
2216

2317
class HouLiFilterCached : HouLiFilter {

‎cpp/lib/Naive.cpp‎

Lines changed: 69 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -45,9 +45,9 @@ void Naive::init(std::string_view equilibriumName) {
4545

4646
// aPar, uekPar, ne
4747
g.for_each_kxky([&](Dim kx, Dim ky) {
48-
moments_K(kx, ky, N_E) = nonlinear::phiInv(phi_K(kx, ky), kPerp2(kx, ky));
48+
moments_K(kx, ky, N_E) = nonlinear::phiInv(phi_K(kx, ky), g.kPerp2(kx, ky));
4949
moments_K(kx, ky, A_PAR) = aParEq_K(kx, ky);
50-
ueKPar_K(kx, ky) = -kPerp2(kx, ky) * moments_K(kx, ky, A_PAR);
50+
ueKPar_K(kx, ky) = -g.kPerp2(kx, ky) * moments_K(kx, ky, A_PAR);
5151
});
5252
br.derivatives(phi_K, dPhi);
5353
br.derivatives(ueKPar_K, dUEKPar);
@@ -121,7 +121,7 @@ void Naive::run(Dim N, Dim saveInterval) {
121121
dt / 2.0 * (1 + exp_nu(kx, ky, hyper.nu_2, dt)) * GM_Nonlinear_K(kx, ky, N_E);
122122

123123
GM_Nonlinear_K(kx, ky, A_PAR) =
124-
nonlinear::A(bracketAParPhiG2Ne_K(kx, ky), bracketUEParPhi_K(kx, ky), kPerp2(kx, ky));
124+
nonlinear::A(bracketAParPhiG2Ne_K(kx, ky), bracketUEParPhi_K(kx, ky), g.kPerp2(kx, ky));
125125
GM_K_Star(kx, ky, A_PAR) =
126126
exp_eta(kx, ky, hyper.eta2, dt) * moments_K(kx, ky, A_PAR) +
127127
dt / 2.0 * (1 + exp_eta(kx, ky, hyper.eta2, dt)) * GM_Nonlinear_K(kx, ky, A_PAR) +
@@ -138,7 +138,7 @@ void Naive::run(Dim N, Dim saveInterval) {
138138
auto bracketPhiGLast_K = br.halfBracket(dPhi, Grid::sliceXY(dGM, LAST));
139139
auto bracketAParGLast_K = br.halfBracket(Grid::sliceXY(dGM, A_PAR), Grid::sliceXY(dGM, LAST));
140140
g.for_each_kxky([&](Dim kx, Dim ky) {
141-
bracketAParGLast_K(kx, ky) *= nonlinear::GLastBracketFactor(g.M, kPerp2(kx, ky), hyper);
141+
bracketAParGLast_K(kx, ky) *= nonlinear::GLastBracketFactor(g.M, g.kPerp2(kx, ky), hyper);
142142
bracketAParGLast_K(kx, ky) += rhoS / de * std::sqrt(LAST) * moments_K(kx, ky, LAST - 1);
143143
// TODO Viriato adds this after the derivative
144144
});
@@ -192,8 +192,8 @@ void Naive::run(Dim N, Dim saveInterval) {
192192
g.for_each_kxky([&](Dim kx, Dim ky) {
193193
// set to 0 for (kx, ky)=(0,0)
194194
phi_K_New(kx, ky) =
195-
((kx | ky) == 0) ? 0 : nonlinear::phi(GM_K_Star(kx, ky, N_E), kPerp2(kx, ky));
196-
ueKPar_K_New(kx, ky) = -kPerp2(kx, ky) * GM_K_Star(kx, ky, A_PAR);
195+
((kx | ky) == 0) ? 0 : nonlinear::phi(GM_K_Star(kx, ky, N_E), g.kPerp2(kx, ky));
196+
ueKPar_K_New(kx, ky) = -g.kPerp2(kx, ky) * GM_K_Star(kx, ky, A_PAR);
197197
});
198198

199199
auto dPhi_Loop = g.dBufXY(), dUEKPar_Loop = g.dBufXY();
@@ -213,7 +213,7 @@ void Naive::run(Dim N, Dim saveInterval) {
213213
auto guessAPar_K = g.cBufXY(), semiImplicitOperator = g.cBufXY();
214214
g.for_each_kxky([&](Dim kx, Dim ky) {
215215
guessAPar_K(kx, ky) = moments_K(kx, ky, A_PAR);
216-
semiImplicitOperator(kx, ky) = nonlinear::semiImplicitOp(dt, bPerpMax, aa0, kPerp2(kx, ky));
216+
semiImplicitOperator(kx, ky) = nonlinear::semiImplicitOp(dt, bPerpMax, aa0, g.kPerp2(kx, ky));
217217
});
218218
semiImplicitOperator(0, 0) = 0;
219219

@@ -249,7 +249,7 @@ void Naive::run(Dim N, Dim saveInterval) {
249249
Real sumAParRelError = 0;
250250
g.for_each_kxky([&](Dim kx, Dim ky) {
251251
GM_Nonlinear_K_Loop(kx, ky, A_PAR) = nonlinear::A(
252-
bracketAParPhiG2Ne_K_Loop(kx, ky), bracketUEParPhi_K_Loop(kx, ky), kPerp2(kx, ky));
252+
bracketAParPhiG2Ne_K_Loop(kx, ky), bracketUEParPhi_K_Loop(kx, ky), g.kPerp2(kx, ky));
253253
// TODO(OPT) reuse star
254254
momentsNew_K(kx, ky, A_PAR) =
255255
1.0 / (1.0 + semiImplicitOperator(kx, ky) / 4.0) *
@@ -258,7 +258,7 @@ void Naive::run(Dim N, Dim saveInterval) {
258258
dt / 2.0 * GM_Nonlinear_K_Loop(kx, ky, A_PAR) +
259259
(1.0 - exp_eta(kx, ky, hyper.eta2, dt)) * aParEq_K(kx, ky) +
260260
semiImplicitOperator(kx, ky) / 4.0 * guessAPar_K(kx, ky));
261-
ueKPar_K_New(kx, ky) = -kPerp2(kx, ky) * momentsNew_K(kx, ky, A_PAR);
261+
ueKPar_K_New(kx, ky) = -g.kPerp2(kx, ky) * momentsNew_K(kx, ky, A_PAR);
262262

263263
sumAParRelError += std::norm(momentsNew_K(kx, ky, A_PAR) - moments_K(kx, ky, A_PAR));
264264
});
@@ -291,7 +291,7 @@ void Naive::run(Dim N, Dim saveInterval) {
291291
dt / 2.0 * GM_Nonlinear_K_Loop(kx, ky, N_E);
292292

293293
phi_K_New(kx, ky) =
294-
(kx | ky) == 0 ? 0 : nonlinear::phi(momentsNew_K(kx, ky, N_E), kPerp2(kx, ky));
294+
(kx | ky) == 0 ? 0 : nonlinear::phi(momentsNew_K(kx, ky, N_E), g.kPerp2(kx, ky));
295295
});
296296

297297
br.derivatives(phi_K_New, dPhi_Loop);
@@ -347,7 +347,7 @@ void Naive::run(Dim N, Dim saveInterval) {
347347
br.halfBracket(Grid::sliceXY(dGM_Loop, A_PAR), Grid::sliceXY(dGM_Loop, LAST));
348348
g.for_each_kxky([&](Dim kx, Dim ky) {
349349
bracketAParGLast_K_Loop(kx, ky) *=
350-
nonlinear::GLastBracketFactor(g.M, kPerp2(kx, ky), hyper);
350+
nonlinear::GLastBracketFactor(g.M, g.kPerp2(kx, ky), hyper);
351351
bracketAParGLast_K_Loop(kx, ky) +=
352352
rhoS / de * std::sqrt(LAST) * momentsNew_K(kx, ky, LAST - 1);
353353
// Note: Viriato adds this after derivative, but can be distributed
@@ -479,16 +479,70 @@ Naive::Buf::R_XY Naive::getMoment(Dim m) const {
479479
return out;
480480
}
481481

482+
Real Naive::getTimestep(DxDy<View::R_XY> dPhi, DxDy<View::R_XY> dNE, DxDy<View::R_XY> dAPar) {
483+
// compute flows
484+
DxDy<View::R_XY> ve, b;
485+
Real vyMax{0}, vxMax{0}, bxMax{0}, byMax{0};
486+
bPerpMax = 0;
487+
488+
// Note that this is minus in Viriato, but we don't care because we're taking the absolute value
489+
// anyway.
490+
ve.DX = dPhi.DY;
491+
ve.DY = dPhi.DX;
492+
b.DX = dAPar.DY;
493+
b.DY = dAPar.DX;
494+
495+
g.for_each_xy([&](Dim x, Dim y) {
496+
bxMax = std::max(bxMax, std::abs(b.DX(x, y)));
497+
byMax = std::max(byMax, std::abs(b.DY(x, y)));
498+
bPerpMax = std::max(bPerpMax, std::sqrt(b.DX(x, y) * b.DX(x, y) + b.DY(x, y) * b.DY(x, y)));
499+
vxMax = std::max(vxMax, std::abs(ve.DX(x, y)));
500+
vyMax = std::max(vyMax, std::abs(ve.DY(x, y)));
501+
if (rhoI >= smallRhoI) {
502+
vxMax = std::max(vxMax, rhoS * rhoS * std::abs(dNE.DX(x, y)));
503+
vyMax = std::max(vyMax, rhoS * rhoS * std::abs(dNE.DY(x, y)));
504+
}
505+
});
506+
507+
Real kperpDum2 = std::pow(g.ky_(g.KY / 2), 2) + std::pow(Real(g.KX), 2);
508+
Real omegaKaw;
509+
if (rhoI < smallRhoI) {
510+
omegaKaw = std::sqrt(1.0 + kperpDum2 * (3.0 / 4.0 * rhoI * rhoI + rhoS * rhoS)) *
511+
g.ky_(g.KY / 2) * bPerpMax / (1.0 + kperpDum2 * de * de);
512+
} else {
513+
omegaKaw =
514+
std::sqrt(kperpDum2 *
515+
(rhoS * rhoS - rhoI * rhoI / (Gamma0(0.5 * kperpDum2 * rhoI * rhoI) - 1.0))) *
516+
g.ky_(g.KY / 2 + 1) * bPerpMax / std::sqrt(1.0 + kperpDum2 * de * de);
517+
}
518+
519+
Real dx = lx / Real(g.X), dy = ly / Real(g.Y);
520+
521+
Real CFLFlow;
522+
if (g.M > 2) {
523+
CFLFlow = std::min({dx / vxMax, dy / vyMax, 2.0 / omegaKaw,
524+
std::min(dx / bxMax, dy / byMax) / (rhoS / de) / std::sqrt(LAST)});
525+
} else {
526+
CFLFlow = std::min({dx / vxMax, dy / vyMax, 2.0 / omegaKaw, dx / bxMax, dy / byMax});
527+
}
528+
529+
spdlog::debug("vxmax: {}, vymax: {}, bxmax: {}, bymax: {}", vxMax, vyMax, bxMax, byMax);
530+
spdlog::debug("bperp_max: {}, omegakaw: {}, CFLFlow: {}", bPerpMax, omegaKaw, CFLFlow);
531+
spdlog::debug("Calculated timestep: {}", CFLFrac * CFLFlow);
532+
533+
return CFLFrac * CFLFlow;
534+
}
535+
482536
Naive::Energies Naive::calculateEnergies() const {
483537
Energies e{};
484538
g.for_each_kxky([&](Dim kx, Dim ky) {
485539
Real const factor = kx == 0 ? 0.5 : 1.0;
486-
e.magnetic += factor * kPerp2(kx, ky) * std::norm(moments_K(kx, ky, A_PAR));
540+
e.magnetic += factor * g.kPerp2(kx, ky) * std::norm(moments_K(kx, ky, A_PAR));
487541
if (rhoI < smallRhoI) {
488-
e.kinetic += factor * kPerp2(kx, ky) * std::norm(phi_K(kx, ky));
542+
e.kinetic += factor * g.kPerp2(kx, ky) * std::norm(phi_K(kx, ky));
489543
} else {
490-
e.kinetic -= factor * 1.0 / (rhoI * rhoI) * (Gamma0(kPerp2(kx, ky) * rhoI * rhoI / 2.0) - 1) *
491-
std::norm(phi_K(kx, ky));
544+
e.kinetic -= factor * 1.0 / (rhoI * rhoI) *
545+
(Gamma0(g.kPerp2(kx, ky) * rhoI * rhoI / 2.0) - 1) * std::norm(phi_K(kx, ky));
492546
}
493547
});
494548

‎cpp/lib/Naive.hpp‎

Lines changed: 5 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include "constants.hpp"
99
#include "debug.hpp"
1010
#include "grid.hpp"
11+
#include "hyper.hpp"
1112
#include "nonlinears.hpp"
1213

1314
#include <fftw-cpp/fftw-cpp.h>
@@ -96,87 +97,23 @@ class Naive : public ahr::HermiteRunner {
9697
// TODO other file/class
9798
// =================
9899

99-
[[nodiscard]] Real ky_(Dim ky) const {
100-
return (ky <= (g.KY / 2) ? Real(ky) : Real(ky) - Real(g.KY)) * Real(lx) / Real(ly);
101-
}
102-
[[nodiscard]] Real kx_(Dim kx) const { return Real(kx); }
103-
104-
[[nodiscard]] Real kPerp2(Dim kx, Dim ky) const {
105-
auto dkx = kx_(kx), dky = ky_(ky);
106-
return dkx * dkx + dky * dky;
107-
}
108-
109-
[[nodiscard]] Real kPerp(Dim kx, Dim ky) const { return std::sqrt(kPerp2(kx, ky)); }
110-
111100
[[nodiscard]] Real exp_nu(Dim kx, Dim ky, Real nu2, Real dt) const {
112-
return std::exp(-(nu * kPerp2(kx, ky) + nu2 * std::pow(kPerp2(kx, ky), hyper_order)) * dt);
101+
return std::exp(-(nu * g.kPerp2(kx, ky) + nu2 * std::pow(g.kPerp2(kx, ky), hyper_order)) * dt);
113102
}
114103

115104
[[nodiscard]] Real exp_gm(Dim m, Real hyper_nuei, Real dt) const {
116105
return exp(-(Real(m) * nu_ei + std::pow(m, 2 * hyper_morder) * hyper_nuei) * dt);
117106
}
118107

119108
[[nodiscard]] Real exp_eta(Dim kx, Dim ky, Real res2, Real dt) const {
120-
return std::exp(-(res * kPerp2(kx, ky) + res2 * std::pow(kPerp2(kx, ky), hyper_order)) * dt /
121-
(1.0 + kPerp2(kx, ky) * de * de));
109+
return std::exp(-(res * g.kPerp2(kx, ky) + res2 * std::pow(g.kPerp2(kx, ky), hyper_order)) *
110+
dt / (1.0 + g.kPerp2(kx, ky) * de * de));
122111
}
123112

124113
/// getTimestep calculates flows and magnetic fields to determine a dt.
125114
/// It also updates bPerpMax in the process.
126115
[[nodiscard]] Real getTimestep(DxDy<View::R_XY> dPhi, DxDy<View::R_XY> dNE,
127-
DxDy<View::R_XY> dAPar) {
128-
// compute flows
129-
DxDy<View::R_XY> ve, b;
130-
Real vyMax{0}, vxMax{0}, bxMax{0}, byMax{0};
131-
bPerpMax = 0;
132-
133-
// Note that this is minus in Viriato, but we don't care because we're taking the absolute value
134-
// anyway.
135-
ve.DX = dPhi.DY;
136-
ve.DY = dPhi.DX;
137-
b.DX = dAPar.DY;
138-
b.DY = dAPar.DX;
139-
140-
g.for_each_xy([&](Dim x, Dim y) {
141-
bxMax = std::max(bxMax, std::abs(b.DX(x, y)));
142-
byMax = std::max(byMax, std::abs(b.DY(x, y)));
143-
bPerpMax = std::max(bPerpMax, std::sqrt(b.DX(x, y) * b.DX(x, y) + b.DY(x, y) * b.DY(x, y)));
144-
vxMax = std::max(vxMax, std::abs(ve.DX(x, y)));
145-
vyMax = std::max(vyMax, std::abs(ve.DY(x, y)));
146-
if (rhoI >= smallRhoI) {
147-
vxMax = std::max(vxMax, rhoS * rhoS * std::abs(dNE.DX(x, y)));
148-
vyMax = std::max(vyMax, rhoS * rhoS * std::abs(dNE.DY(x, y)));
149-
}
150-
});
151-
152-
Real kperpDum2 = std::pow(ky_(g.KY / 2), 2) + std::pow(Real(g.KX), 2);
153-
Real omegaKaw;
154-
if (rhoI < smallRhoI) {
155-
omegaKaw = std::sqrt(1.0 + kperpDum2 * (3.0 / 4.0 * rhoI * rhoI + rhoS * rhoS)) *
156-
ky_(g.KY / 2) * bPerpMax / (1.0 + kperpDum2 * de * de);
157-
} else {
158-
omegaKaw =
159-
std::sqrt(kperpDum2 *
160-
(rhoS * rhoS - rhoI * rhoI / (Gamma0(0.5 * kperpDum2 * rhoI * rhoI) - 1.0))) *
161-
ky_(g.KY / 2 + 1) * bPerpMax / std::sqrt(1.0 + kperpDum2 * de * de);
162-
}
163-
164-
Real dx = lx / Real(g.X), dy = ly / Real(g.Y);
165-
166-
Real CFLFlow;
167-
if (g.M > 2) {
168-
CFLFlow = std::min({dx / vxMax, dy / vyMax, 2.0 / omegaKaw,
169-
std::min(dx / bxMax, dy / byMax) / (rhoS / de) / std::sqrt(LAST)});
170-
} else {
171-
CFLFlow = std::min({dx / vxMax, dy / vyMax, 2.0 / omegaKaw, dx / bxMax, dy / byMax});
172-
}
173-
174-
spdlog::debug("vxmax: {}, vymax: {}, bxmax: {}, bymax: {}", vxMax, vyMax, bxMax, byMax);
175-
spdlog::debug("bperp_max: {}, omegakaw: {}, CFLFlow: {}", bPerpMax, omegaKaw, CFLFlow);
176-
spdlog::debug("Calculated timestep: {}", CFLFrac * CFLFlow);
177-
178-
return CFLFrac * CFLFlow;
179-
}
116+
DxDy<View::R_XY> dAPar);
180117

181118
public:
182119
struct Energies {

‎cpp/lib/constants.hpp‎

Lines changed: 0 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,5 @@
11
#pragma once
22

3-
#include "grid.hpp"
43
#include "typedefs.hpp"
54

65
namespace ahr {
@@ -44,26 +43,4 @@ constexpr Dim hyper_morder = 3;
4443

4544
/// equil
4645
constexpr Real a0 = 1.0;
47-
48-
struct HyperCoefficients {
49-
Real nu_g, nu_2, eta2, nu_ei;
50-
51-
static HyperCoefficients calculate(Real dt, Grid const &g) {
52-
Real kPerpMax2 = std::pow(g.KX, 2) + std::pow(Real(g.KY) / 2, 2);
53-
54-
HyperCoefficients ret{};
55-
ret.nu_g = hyper_coef_g / dt / std::pow(kPerpMax2, hyper_order_g);
56-
ret.nu_2 = hyper_coef / dt / std::pow(kPerpMax2, hyper_order);
57-
if (kPerpMax2 * de * de > 1) {
58-
ret.eta2 = hyper_coef / dt / std::pow(kPerpMax2, hyper_order - 1) * de * de;
59-
} else {
60-
ret.eta2 = hyper_coef / dt / std::pow(kPerpMax2, hyper_order);
61-
}
62-
63-
ret.nu_ei = hyperm_coef / dt / std::pow(g.M, 2 * hyper_morder);
64-
65-
return ret;
66-
}
67-
};
68-
6946
} // namespace ahr

‎cpp/lib/grid.hpp‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#pragma once
22

3+
#include "constants.hpp"
34
#include "typedefs.hpp"
45

56
#include <fftw-cpp/fftw-cpp.h>
@@ -117,6 +118,18 @@ struct Grid {
117118
}
118119
}
119120
}
121+
122+
[[nodiscard]] Real ky_(Dim ky) const {
123+
return (ky <= (KY / 2) ? Real(ky) : Real(ky) - Real(KY)) * Real(lx) / Real(ly);
124+
}
125+
[[nodiscard]] Real kx_(Dim kx) const { return Real(kx); }
126+
127+
[[nodiscard]] Real kPerp2(Dim kx, Dim ky) const {
128+
auto dkx = kx_(kx), dky = ky_(ky);
129+
return dkx * dkx + dky * dky;
130+
}
131+
132+
[[nodiscard]] Real kPerp(Dim kx, Dim ky) const { return std::sqrt(kPerp2(kx, ky)); }
120133
};
121134

122135
// Options for various buffer types:

‎cpp/lib/hyper.hpp‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
#pragma once
2+
#include "constants.hpp"
3+
#include "grid.hpp"
4+
5+
namespace ahr {
6+
struct HyperCoefficients {
7+
Real nu_g, nu_2, eta2, nu_ei;
8+
9+
static HyperCoefficients calculate(Real dt, Grid const &g) {
10+
Real kPerpMax2 = std::pow(g.KX, 2) + std::pow(Real(g.KY) / 2, 2);
11+
12+
HyperCoefficients ret{};
13+
ret.nu_g = hyper_coef_g / dt / std::pow(kPerpMax2, hyper_order_g);
14+
ret.nu_2 = hyper_coef / dt / std::pow(kPerpMax2, hyper_order);
15+
if (kPerpMax2 * de * de > 1) {
16+
ret.eta2 = hyper_coef / dt / std::pow(kPerpMax2, hyper_order - 1) * de * de;
17+
} else {
18+
ret.eta2 = hyper_coef / dt / std::pow(kPerpMax2, hyper_order);
19+
}
20+
21+
ret.nu_ei = hyperm_coef / dt / std::pow(g.M, 2 * hyper_morder);
22+
23+
return ret;
24+
}
25+
};
26+
} // namespace ahr

0 commit comments

Comments
 (0)