Skip to content

Commit 222f4ef

Browse files
UCaromelnicolasaunai
UCaromel
authored and
nicolasaunai
committed
weno-z
1 parent 6dd8c6f commit 222f4ef

File tree

4 files changed

+135
-11
lines changed

4 files changed

+135
-11
lines changed

‎src/core/numerics/reconstructions/wenoz.hpp

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,113 @@
1+
#ifndef CORE_NUMERICS_RECONSTRUCTION_WENOZ_HPP
2+
#define CORE_NUMERICS_RECONSTRUCTION_WENOZ_HPP
3+
4+
#include "core/data/vecfield/vecfield_component.hpp"
5+
#include "core/utilities/index/index.hpp"
6+
#include <utility>
7+
8+
namespace PHARE::core
9+
{
10+
template<typename GridLayout, typename SlopeLimiter = void>
11+
class WENOZReconstruction
12+
{
13+
public:
14+
using GridLayout_t = GridLayout;
15+
16+
template<auto direction, typename Field>
17+
static auto reconstruct(Field const& F, MeshIndex<Field::dimension> index)
18+
{
19+
auto u_3
20+
= F(GridLayout::template previous<direction>(GridLayout::template previous<direction>(
21+
GridLayout::template previous<direction>(index))));
22+
auto u_2 = F(GridLayout::template previous<direction>(
23+
GridLayout::template previous<direction>(index)));
24+
auto u_1 = F(GridLayout::template previous<direction>(index));
25+
auto u = F(index);
26+
auto u1 = F(GridLayout::template next<direction>(index));
27+
auto u2
28+
= F(GridLayout::template next<direction>(GridLayout::template next<direction>(index)));
29+
30+
return std::make_pair(recons_wenoz_L_(u_3, u_2, u_1, u, u1),
31+
recons_wenoz_R_(u_2, u_1, u, u1, u2));
32+
}
33+
34+
template<auto direction, typename Field>
35+
static auto center_reconstruct(Field const& U, MeshIndex<Field::dimension> index,
36+
auto projection)
37+
{
38+
auto u_3 = GridLayout::project(
39+
U,
40+
GridLayout::template previous<direction>(GridLayout::template previous<direction>(
41+
GridLayout::template previous<direction>(index))),
42+
projection);
43+
auto u_2 = GridLayout::project(U,
44+
GridLayout::template previous<direction>(
45+
GridLayout::template previous<direction>(index)),
46+
projection);
47+
auto u_1
48+
= GridLayout::project(U, GridLayout::template previous<direction>(index), projection);
49+
auto u = GridLayout::project(U, index, projection);
50+
auto u1 = GridLayout::project(U, GridLayout::template next<direction>(index), projection);
51+
auto u2 = GridLayout::project(
52+
U, GridLayout::template next<direction>(GridLayout::template next<direction>(index)),
53+
projection);
54+
55+
return std::make_pair(recons_wenoz_L_(u_3, u_2, u_1, u, u1),
56+
recons_wenoz_R_(u_2, u_1, u, u1, u2));
57+
}
58+
59+
private:
60+
static auto recons_wenoz_L_(auto const ull, auto const ul, auto const u, auto const ur,
61+
auto const urr)
62+
{
63+
static constexpr auto dL0 = 1. / 10.;
64+
static constexpr auto dL1 = 3. / 5.;
65+
static constexpr auto dL2 = 3. / 10.;
66+
67+
auto const [wL0, wL1, wL2] = compute_wenoz_weights(ull, ul, u, ur, urr, dL0, dL1, dL2);
68+
69+
return wL0 * ((1. / 3.) * ull - (7. / 6.) * ul + (11. / 6.) * u)
70+
+ wL1 * (-(1. / 6.) * ul + (5. / 6.) * u + (1. / 3.) * ur)
71+
+ wL2 * ((1. / 3.) * u + (5. / 6.) * ur - (1. / 6.) * urr);
72+
}
73+
74+
static auto recons_wenoz_R_(auto const ull, auto const ul, auto const u, auto const ur,
75+
auto const urr)
76+
{
77+
static constexpr auto dR0 = 3. / 10.;
78+
static constexpr auto dR1 = 3. / 5.;
79+
static constexpr auto dR2 = 1. / 10.;
80+
81+
auto const [wR0, wR1, wR2] = compute_wenoz_weights(ull, ul, u, ur, urr, dR0, dR1, dR2);
82+
83+
return wR0 * ((1. / 3.) * u + (5. / 6.) * ul - (1. / 6.) * ull)
84+
+ wR1 * (-(1. / 6.) * ur + (5. / 6.) * u + (1. / 3.) * ul)
85+
+ wR2 * ((1. / 3.) * urr - (7. / 6.) * ur + (11. / 6.) * u);
86+
}
87+
88+
static auto compute_wenoz_weights(auto const ull, auto const ul, auto const u, auto const ur,
89+
auto const urr, auto const d0, auto const d1, auto const d2)
90+
{
91+
constexpr auto eps = 1.e-40;
92+
93+
auto const beta0 = (13. / 12.) * (ull - 2. * ul + u) * (ull - 2. * ul + u)
94+
+ (1. / 4.) * (ull - 4. * ul + 3. * u) * (ull - 4. * ul + 3. * u);
95+
auto const beta1 = (13. / 12.) * (ul - 2. * u + ur) * (ul - 2. * u + ur)
96+
+ (1. / 4.) * (ul - ur) * (ul - ur);
97+
auto const beta2 = (13. / 12.) * (u - 2. * ur + urr) * (u - 2. * ur + urr)
98+
+ (1. / 4.) * (3. * u - 4. * ur + urr) * (3. * u - 4. * ur + urr);
99+
100+
auto const tau5 = std::abs(beta0 - beta2);
101+
102+
auto const alpha0 = d0 * (1. + tau5 / (beta0 + eps));
103+
auto const alpha1 = d1 * (1. + tau5 / (beta1 + eps));
104+
auto const alpha2 = d2 * (1. + tau5 / (beta2 + eps));
105+
106+
auto const sum_alpha = alpha0 + alpha1 + alpha2;
107+
return std::make_tuple(alpha0 / sum_alpha, alpha1 / sum_alpha, alpha2 / sum_alpha);
108+
}
109+
};
110+
111+
} // namespace PHARE::core
112+
113+
#endif

‎src/python3/pyMHD.cpp

Lines changed: 17 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include <pybind11/pybind11.h>
66
#include <pybind11/stl.h>
77

8+
#include "core/numerics/reconstructions/wenoz.hpp"
89
#include "core/utilities/types.hpp"
910
#include "tests/core/numerics/mock_mhd_simulator/test_mhd_solver.hpp"
1011

@@ -57,7 +58,7 @@ template<std::size_t Constant>
5758
using InterpConst = PHARE::core::InterpConst<Constant>;
5859

5960
enum class TimeIntegratorType { Euler, TVDRK2, TVDRK3, count };
60-
enum class ReconstructionType { Constant, Linear, WENO3, count };
61+
enum class ReconstructionType { Constant, Linear, WENO3, WENOZ, count };
6162
enum class SlopeLimiterType { VanLeer, MinMod, count };
6263
enum class RiemannSolverType { Rusanov, HLL, count };
6364

@@ -115,6 +116,13 @@ struct ReconstructionSelector<ReconstructionType::WENO3>
115116
using type = WENO3Reconstruction<GridLayout, SlopeLimiter>;
116117
};
117118

119+
template<>
120+
struct ReconstructionSelector<ReconstructionType::WENOZ>
121+
{
122+
template<typename GridLayout, typename SlopeLimiter>
123+
using type = WENOZReconstruction<GridLayout, SlopeLimiter>;
124+
};
125+
118126
template<ReconstructionType R, SlopeLimiterType S>
119127
struct SlopeLimiterSelector
120128
{
@@ -242,7 +250,8 @@ void registerSimulatorVariants(py::module& m)
242250
+ std::string("_")
243251
+ (rc == ReconstructionType::Constant ? "constant"
244252
: rc == ReconstructionType::Linear ? "linear"
245-
: "weno3")
253+
: rc == ReconstructionType::WENO3 ? "weno3"
254+
: "wenoz")
246255
+ std::string("_")
247256
+ (sl == SlopeLimiterType::VanLeer ? "vanleer" : "minmod")
248257

@@ -263,7 +272,9 @@ void registerSimulatorVariants(py::module& m)
263272
: ti == TimeIntegratorType::TVDRK2 ? "tvdrk2"
264273
: "tvdrk3")
265274
+ std::string("_")
266-
+ (rc == ReconstructionType::Constant ? "constant" : "weno3")
275+
+ (rc == ReconstructionType::Constant ? "constant"
276+
: rc == ReconstructionType::WENO3 ? "weno3"
277+
: "wenoz")
267278
+ std::string("_")
268279
+ (rs == RiemannSolverType::Rusanov ? "rusanov" : "hll")
269280
+ (hall ? "_hall" : "") + (res ? "_res" : "")
@@ -296,7 +307,7 @@ PYBIND11_MODULE(pyMHD, m)
296307
/* SlopeLimiterType::count, RiemannSolverType::Rusanov, false, false,*/
297308
/* false>::registerVariant("2_1_tvdrk2_constant_rusanov", m);*/
298309

299-
Registerer<2, 1, TimeIntegratorType::TVDRK3, ReconstructionType::WENO3, SlopeLimiterType::count,
300-
RiemannSolverType::Rusanov, true, false,
301-
false>::registerVariant("2_1_tvdrk3_weno3_rusanov_hall", m);
310+
Registerer<2, 2, TimeIntegratorType::TVDRK3, ReconstructionType::WENOZ, SlopeLimiterType::count,
311+
RiemannSolverType::Rusanov, false, false,
312+
false>::registerVariant("2_2_tvdrk3_wenoz_rusanov", m);
302313
}

‎tests/functional/mhd_harris/harris.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ def config():
1010

1111
sim = s.Simulation(
1212
ndim=2,
13-
order=1,
13+
order=2,
1414
timestep=0.002,
1515
final_time=10,
1616
cells=cells,
@@ -19,7 +19,7 @@ def config():
1919
eta=0.0,
2020
nu=0.0,
2121
gamma=5.0 / 3.0,
22-
reconstruction="weno3",
22+
reconstruction="wenoz",
2323
limiter="",
2424
riemann="rusanov",
2525
time_integrator="tvdrk3",

‎tests/functional/mhd_orszagtang/orszag_tang.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,16 +7,16 @@
77
def config():
88
sim = s.Simulation(
99
ndim=2,
10-
order=1,
10+
order=2,
1111
timestep=0.0014,
12-
final_time=0.5,
12+
final_time=3,
1313
cells=(128, 128),
1414
dl=(1.0 / 128.0, 1.0 / 128.0),
1515
origin=(0.0, 0.0),
1616
eta=0.0,
1717
nu=0.0,
1818
gamma=5.0 / 3.0,
19-
reconstruction="weno3",
19+
reconstruction="wenoz",
2020
limiter="",
2121
riemann="rusanov",
2222
time_integrator="tvdrk3",

0 commit comments

Comments
 (0)