DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
SpecialFields.hpp
Go to the documentation of this file.
1/** @file SpecialFields.hpp
2 * @brief Analytic isentropic-vortex solutions for inviscid accuracy verification.
3 *
4 * Provides three variants of the classical isentropic vortex used as exact
5 * solutions for the Euler equations:
6 * - IsentropicVortex10 — moving vortex on a [0,10]^2 periodic domain.
7 * - IsentropicVortex30 — moving vortex on a [-10,20]^2 periodic domain (period 30).
8 * - IsentropicVortexCent — stationary vortex centered at the origin (no periodicity).
9 *
10 * All three share the same isentropic vortex formula:
11 * - Temperature perturbation: dT = -(gamma-1)/(8*gamma*pi^2) * chi^2 * exp(1-r^2)
12 * - Velocity perturbation (solid-body rotation):
13 * - du_x = -chi/(2*pi) * exp((1-r^2)/2) * (y - y_c)
14 * - du_y = +chi/(2*pi) * exp((1-r^2)/2) * (x - x_c)
15 * - T = 1 + dT, rho = T^(1/(gamma-1)), p = T * rho
16 *
17 * The returned state vector is in conservative variables: (rho, rho*u, rho*v, [rho*w,] E).
18 */
19#pragma once
20
21#include "EulerEvaluator.hpp"
22
23namespace DNDS::Euler::SpecialFields
24{
25 /**
26 * @brief Analytic isentropic vortex on a [0,10] x [0,10] periodic domain.
27 *
28 * Background flow (u,v) = (1,1), vortex center at (5,5). The domain wraps
29 * with period 10 in both x and y via float_mod, so the vortex convects
30 * diagonally and re-enters from the opposite corner. Used for inviscid
31 * accuracy verification (spatial and temporal order-of-accuracy studies).
32 *
33 * @tparam model EulerModel selecting the equation set (default NS = 2-D Navier-Stokes).
34 * @param eval EulerEvaluator providing gas properties (gamma).
35 * @param x Physical-space coordinates at which to evaluate the solution.
36 * @param t Current physical time (the vortex translates at speed (1,1)*t).
37 * @param cnVars Number of conservative variables in the returned vector.
38 * @param chi Vortex strength parameter (standard test value = 5).
39 * @return Conservative state vector (rho, rho*u_x, rho*u_y, ..., E) of size cnVars.
40 */
41 template <EulerModel model = NS>
42 auto IsentropicVortex10(
43 EulerEvaluator<model> &eval,
44 const Geom::tPoint &x,
45 real t, int cnVars,
46 real chi)
47 {
48 typename EulerEvaluator<model>::TU ret;
49 ret.resize(cnVars);
50
51 real xyc = 5;
52 real gamma = eval.settings.idealGasProperty.gamma;
53 Geom::tPoint pPhysics = x;
54 pPhysics[0] = float_mod(pPhysics[0] - t, 10);
55 pPhysics[1] = float_mod(pPhysics[1] - t, 10);
56 real r = std::sqrt(sqr(pPhysics(0) - xyc) + sqr(pPhysics(1) - xyc));
57 real dT = -(gamma - 1) / (8 * gamma * sqr(pi)) * sqr(chi) * std::exp(1 - sqr(r));
58 real dux = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * -(pPhysics(1) - xyc);
59 real duy = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * +(pPhysics(0) - xyc);
60 real T = dT + 1;
61 real ux = dux + 1;
62 real uy = duy + 1;
63 real S = 1;
64 real rho = std::pow(T / S, 1 / (gamma - 1));
65 real p = T * rho;
66
67 real E = 0.5 * (sqr(ux) + sqr(uy)) * rho + p / (gamma - 1);
68
69 ret.setZero();
70 ret(0) = rho;
71 ret(1) = rho * ux;
72 ret(2) = rho * uy;
73 ret(EulerEvaluator<model>::dim + 1) = E;
74 return ret;
75 }
76
77 /**
78 * @brief Analytic isentropic vortex on a [-10,20] x [-10,20] periodic domain (period 30).
79 *
80 * Same vortex formula as IsentropicVortex10 but on a larger domain with
81 * period 30. Background flow (u,v) = (1,1), vortex center at (5,5),
82 * chi = 5 (hardcoded). Useful when more space around the vortex core is
83 * needed to reduce periodic-image interactions in convergence studies.
84 *
85 * @tparam model EulerModel selecting the equation set (default NS).
86 * @param eval EulerEvaluator providing gas properties (gamma).
87 * @param x Physical-space coordinates at which to evaluate the solution.
88 * @param t Current physical time (vortex translates at speed (1,1)*t).
89 * @param cnVars Number of conservative variables in the returned vector.
90 * @return Conservative state vector (rho, rho*u_x, rho*u_y, ..., E) of size cnVars.
91 */
92 template <EulerModel model = NS>
93 auto IsentropicVortex30(
94 EulerEvaluator<model> &eval,
95 const Geom::tPoint &x,
96 real t, int cnVars)
97 {
98 typename EulerEvaluator<model>::TU ret;
99 ret.resize(cnVars);
100
101 real chi = 5;
102 real xyc = 5;
103 real gamma = eval.settings.idealGasProperty.gamma;
104 Geom::tPoint pPhysics = x;
105 pPhysics[0] = float_mod(pPhysics[0] - t + 10, 30) - 10;
106 pPhysics[1] = float_mod(pPhysics[1] - t + 10, 30) - 10;
107 real r = std::sqrt(sqr(pPhysics(0) - xyc) + sqr(pPhysics(1) - xyc));
108 real dT = -(gamma - 1) / (8 * gamma * sqr(pi)) * sqr(chi) * std::exp(1 - sqr(r));
109 real dux = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * -(pPhysics(1) - xyc);
110 real duy = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * +(pPhysics(0) - xyc);
111 real T = dT + 1;
112 real ux = dux + 1;
113 real uy = duy + 1;
114 real S = 1;
115 real rho = std::pow(T / S, 1 / (gamma - 1));
116 real p = T * rho;
117
118 real E = 0.5 * (sqr(ux) + sqr(uy)) * rho + p / (gamma - 1);
119
120 ret.setZero();
121 ret(0) = rho;
122 ret(1) = rho * ux;
123 ret(2) = rho * uy;
124 ret(EulerEvaluator<model>::dim + 1) = E;
125 return ret;
126 }
127
128 /**
129 * @brief Stationary isentropic vortex centered at the origin with zero background flow.
130 *
131 * Same vortex formula as IsentropicVortex10 but with (u_bg, v_bg) = (0,0) and
132 * the vortex center at the origin (0,0). No periodic wrapping is applied, so
133 * the computational domain should be large enough that the vortex decays
134 * before reaching the boundaries. chi = 5 (hardcoded).
135 *
136 * Used for steady-state vortex preservation tests to verify that the
137 * numerical scheme maintains the stationary vortex without spurious
138 * dissipation or distortion.
139 *
140 * @tparam model EulerModel selecting the equation set (default NS).
141 * @param eval EulerEvaluator providing gas properties (gamma).
142 * @param x Physical-space coordinates at which to evaluate the solution.
143 * @param t Current physical time (unused since the vortex is stationary).
144 * @param cnVars Number of conservative variables in the returned vector.
145 * @return Conservative state vector (rho, rho*u_x, rho*u_y, ..., E) of size cnVars.
146 */
147 template <EulerModel model = NS>
148 auto IsentropicVortexCent(
149 EulerEvaluator<model> &eval,
150 const Geom::tPoint &x,
151 real t, int cnVars)
152 {
153 typename EulerEvaluator<model>::TU ret;
154 ret.resize(cnVars);
155
156 real chi = 5;
157 real xyc = 0; // center is origin
158 real gamma = eval.settings.idealGasProperty.gamma;
159 Geom::tPoint pPhysics = x;
160 // pPhysics[0] = float_mod(pPhysics[0] - t, 10);
161 // pPhysics[1] = float_mod(pPhysics[1] - t, 10);
162 real r = std::sqrt(sqr(pPhysics(0) - xyc) + sqr(pPhysics(1) - xyc));
163 real dT = -(gamma - 1) / (8 * gamma * sqr(pi)) * sqr(chi) * std::exp(1 - sqr(r));
164 real dux = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * -(pPhysics(1) - xyc);
165 real duy = chi / 2 / pi * std::exp((1 - sqr(r)) / 2) * +(pPhysics(0) - xyc);
166 real T = dT + 1;
167 real ux = dux + 0; // no translation
168 real uy = duy + 0;
169 real S = 1;
170 real rho = std::pow(T / S, 1 / (gamma - 1));
171 real p = T * rho;
172
173 real E = 0.5 * (sqr(ux) + sqr(uy)) * rho + p / (gamma - 1);
174
175 ret.setZero();
176 ret(0) = rho;
177 ret(1) = rho * ux;
178 ret(2) = rho * uy;
179 ret(EulerEvaluator<model>::dim + 1) = E;
180 return ret;
181 }
182}
EulerEvaluator.hpp
Core finite-volume evaluator for compressible Navier-Stokes / Euler equations.
DNDS::Euler::EulerEvaluator
Core finite-volume evaluator for compressible Navier-Stokes / Euler equations.
Definition EulerEvaluator.hpp:78
DNDS::Euler::EulerEvaluator::settings
EulerEvaluatorSettings< model > settings
Physics and numerics settings for this evaluator.
Definition EulerEvaluator.hpp:220
DNDS::Euler::SpecialFields::IsentropicVortexCent
auto IsentropicVortexCent(EulerEvaluator< model > &eval, const Geom::tPoint &x, real t, int cnVars)
Stationary isentropic vortex centered at the origin with zero background flow.
Definition SpecialFields.hpp:148
DNDS::Euler::SpecialFields::IsentropicVortex30
auto IsentropicVortex30(EulerEvaluator< model > &eval, const Geom::tPoint &x, real t, int cnVars)
Analytic isentropic vortex on a [-10,20] x [-10,20] periodic domain (period 30).
Definition SpecialFields.hpp:93
DNDS::Euler::SpecialFields::IsentropicVortex10
auto IsentropicVortex10(EulerEvaluator< model > &eval, const Geom::tPoint &x, real t, int cnVars, real chi)
Analytic isentropic vortex on a [0,10] x [0,10] periodic domain.
Definition SpecialFields.hpp:42
DNDS::Geom::tPoint
Eigen::Vector3d tPoint
Definition Geometric.hpp:9
DNDS::pi
DNDS_CONSTANT const real pi
π in double precision (matches DNDS_E_PI macro).
Definition Defines.hpp:199
DNDS::sqr
constexpr T sqr(const T &a)
a * a, constexpr. Works for all arithmetic types.
Definition Defines.hpp:517
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::float_mod
real float_mod(real a, real b)
Floating-point modulo matching Python's % (result has sign of b).
Definition Defines.hpp:580
Eigen::MatrixFMTSafe
Eigen::Matrix wrapper that hides begin/end from fmt.
Definition EigenUtil.hpp:62
r
tVec r(NCells)
x
tVec x(NCells)