Nektar++
TimeRiemann.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TimeRiemann.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:
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include "../RiemannSolvers/RoeSolver.h"
36
38
40
43
44#include "../RiemannSolvers/RoeSolver.h"
45#include "../RiemannSolvers/RoeSolverSIMD.h"
46
47using namespace Nektar;
48
49int main(int argc, char const *argv[])
50{
51
54 LIKWID_MARKER_REGISTER("Riemann");
55 LIKWID_MARKER_REGISTER("RotationTo");
56 LIKWID_MARKER_REGISTER("v_Solve");
57 LIKWID_MARKER_REGISTER("RotationFrom");
58
59 size_t nEle;
60 if (argc < 2)
61 {
62 nEle = 10;
63 }
64 else
65 {
66 nEle = std::stoi(argv[1]);
67 }
68
69 std::cout << "number of faces\t" << nEle
70 << "\t(assuming 4*4 nodes per face)\n";
71
72 // auto riemannSolver = RoeSolver();
73 auto riemannSolver = RoeSolverSIMD();
74 // Setting up parameters for Riemann solver
75 NekDouble gamma = 1.4;
76 riemannSolver.SetParam("gamma",
77 [&gamma]() -> NekDouble & { return gamma; });
78
79 size_t spaceDim = 3;
80 size_t npts = 4 * 4 * nEle; // so that avx spillover loop is engaged
81
82 // Set up locations of velocity vector.
84 vecLocs[0] = Array<OneD, NekDouble>(spaceDim);
85 for (size_t i = 0; i < spaceDim; ++i)
86 {
87 vecLocs[0][i] = 1 + i;
88 }
89 riemannSolver.SetAuxVec(
90 "vecLocs",
91 [&vecLocs]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
92 return vecLocs;
93 });
94
95 // setup normals
96 Array<OneD, Array<OneD, NekDouble>> normals(spaceDim);
97 for (size_t i = 0; i < spaceDim; ++i)
98 {
99 normals[i] = Array<OneD, NekDouble>(npts);
100 }
101 riemannSolver.SetVector(
102 "N", [&normals]() -> const Array<OneD, const Array<OneD, NekDouble>> & {
103 return normals;
104 });
105
106 size_t nFields = spaceDim + 2;
107 Array<OneD, Array<OneD, NekDouble>> fwd(nFields), bwd(nFields),
108 flx(nFields), flxRef(nFields);
109 for (size_t i = 0; i < nFields; ++i)
110 {
111 fwd[i] = Array<OneD, NekDouble>(npts);
112 bwd[i] = Array<OneD, NekDouble>(npts);
113 flx[i] = Array<OneD, NekDouble>(npts);
114 flxRef[i] = Array<OneD, NekDouble>(npts);
115 }
116
117 // up to now it is "boiler plate" code to set up the test
118 // below are the conditions for the test
119
120 for (size_t i = 0; i < npts; ++i)
121 {
122 // density
123 NekDouble rho = 0.9;
124 fwd[0][i] = rho;
125 bwd[0][i] = rho;
126 // x-momentum
127 NekDouble rhou = rho * 1.0;
128 fwd[1][i] = rhou;
129 bwd[1][i] = rhou;
130 // y-momentum
131 NekDouble rhov = rho * 2.0;
132 fwd[2][i] = rhov;
133 bwd[2][i] = rhov;
134 // z-momentum
135 NekDouble rhow = rho * 3.0;
136 fwd[3][i] = rhow;
137 bwd[3][i] = rhow;
138 // energy
139 NekDouble p = 1.0;
140 NekDouble rhoe = p / (gamma - 1.0);
141 NekDouble E =
142 rhoe + 0.5 * (rhou * rhou + rhov * rhov + rhow * rhow) / rho;
143 fwd[nFields - 1][i] = E;
144 bwd[nFields - 1][i] = E;
145 // set face normal along x
146 normals[0][i] = 1.0;
147 // Ref solution
148 flxRef[0][i] = rhou;
149 flxRef[1][i] = rhou * rhou / rho + p;
150 flxRef[2][i] = rhou * rhov / rho;
151 flxRef[3][i] = rhou * rhow / rho;
152 flxRef[nFields - 1][i] = (E + p) * rhou / rho;
153 }
154
155 // number of experiments
156 constexpr size_t experiments = 1 << 12;
157
158 LIKWID_MARKER_START("Riemann");
159 for (size_t j = 0; j < experiments; ++j)
160 {
161 // time
162 riemannSolver.Solve(spaceDim, fwd, bwd, flx);
163 }
164 LIKWID_MARKER_STOP("Riemann");
165 // get likwid events
166 constexpr short CPU_CLK_UNHALTED_REF_id = 2;
167 int nevents{20};
168 std::vector<double> events(nevents);
169 [[maybe_unused]] double time;
170 [[maybe_unused]] int count;
171
172 LIKWID_MARKER_GET("Riemann", &nevents, events.data(), &time, &count);
173 // print out CPE
174 double cpeRiemann = events[CPU_CLK_UNHALTED_REF_id] / npts / experiments;
175 std::cout << "Riemann likwid CPE\t" << cpeRiemann << '\n';
176 //
177 LIKWID_MARKER_GET("RotationTo", &nevents, events.data(), &time, &count);
178 // print out CPE
179 double cpeRotationTo = events[CPU_CLK_UNHALTED_REF_id] / npts / experiments;
180 std::cout << "RotationTo likwid CPE\t" << cpeRotationTo << '\n';
181 //
182 LIKWID_MARKER_GET("v_Solve", &nevents, events.data(), &time, &count);
183 // print out CPE
184 double cpev_Solve = events[CPU_CLK_UNHALTED_REF_id] / npts / experiments;
185 std::cout << "v_Solve likwid CPE\t" << cpev_Solve << '\n';
186 //
187 LIKWID_MARKER_GET("RotationFrom", &nevents, events.data(), &time, &count);
188 // print out CPE
189 double cpeRotationFrom =
190 events[CPU_CLK_UNHALTED_REF_id] / npts / experiments;
191 std::cout << "RotationFrom likwid CPE\t" << cpeRotationFrom << '\n';
192 // avoid opt out
193 std::cout << flx[0][0] << std::endl;
194
196}
#define LIKWID_MARKER_THREADINIT
Definition: Likwid.hpp:43
#define LIKWID_MARKER_START(regionTag)
Definition: Likwid.hpp:46
#define LIKWID_MARKER_CLOSE
Definition: Likwid.hpp:48
#define LIKWID_MARKER_INIT
Definition: Likwid.hpp:42
#define LIKWID_MARKER_REGISTER(regionTag)
Definition: Likwid.hpp:45
#define LIKWID_MARKER_STOP(regionTag)
Definition: Likwid.hpp:47
#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count)
Definition: Likwid.hpp:49
int main(int argc, char const *argv[])
Definition: TimeRiemann.cpp:49
double NekDouble