Peano 4
Loading...
Searching...
No Matches
standaloneRungeKutta.cpph
Go to the documentation of this file.
1#ifndef __KERNELS_STANDALONE_RK_INTEGRATOR_SVEN__
2#define __KERNELS_STANDALONE_RK_INTEGRATOR_SVEN__
3
17#include <algorithm>
18
21#include "computeGradients.cpph"
22
23#include "tarch/la/Vector.h"
24
25namespace kernels {
26namespace aderdg {
27namespace generic {
28namespace c {
29
46template <typename SolverType>
49
50 static constexpr int numberOfVariables = SolverType::NumberOfVariables;
51 static constexpr int numberOfParameters = SolverType::NumberOfParameters;
53 static constexpr int order = SolverType::Order;
54 static constexpr int basisSize = order+1;
55 static constexpr int basisX = basisSize;
56 static constexpr int basisY = basisSize;
57 static constexpr int basisZ = Dimensions==3 ? basisSize : 1;
58 static constexpr int patchSize = basisX*basisY*basisZ*numberOfData;
59
60 // We assume it const here which means that rightHandSide() must has a const qualifier.
61 // You can remove the const if you want, then also remove it in the constructor
62 const SolverType& solver;
63
64 // Cell geometry means the spatial computational domain for this integration
66
67 // Do not mix this up with a CFL number anywhere in ExaHyPE. Here, it is solely for
68 // the RungeKutta timestep determination.
69 const double cfl_number;
70
73
74 void coordinate(double* const pos, const int i) const {
75 pos[i] = offsetOfPatch[i] + sizeOfPatch[i] * (SolverType::nodes[order][i] - 0.5);
76 }
77
78 // Determine the RHS of the equation partial_t Q = S. We used to call this
79 // "FusedSource" in ExaHyPE. Derivatives will be computed here in the DG fashion.
80 void rhs(double* dudtPatch, const double* const QPatch) const {
81 double gradQ[patchSize];
82 kernels::aderdg::generic::c::computeGradQ<SolverType>(gradQ, QPatch, sizeOfPatch);
83
84 const idx5 idx_gradQ(basisX, basisY, basisZ, Dimensions, numberOfVariables);
86
87 double pos[Dimensions];
88
89 for (int iz = 0; iz < basisZ; iz++)
90 for (int iy = 0; iy < basisY; iy++)
91 for (int ix = 0; ix < basisX; ix++) {
92 for(int id=0; id<Dimensions; id++) coordinate(pos, id);
93
94 // I want to pass the position here, so I do not fall back to an ExaHyPE
95 // standard signature such as solver.fusedSource() but use something custom.
96 solver.rightHandSide(
97 pos,
98 QPatch + idx_u(iz,iy,ix,0),
99 gradQ + idx_gradQ(iz,iy,ix,0,0),
100 dudtPatch + idx_u(iz,iy,ix,0)
101 );
102 } // spatial loop
103 } // rhs()
104
105 double largest_eigenvalue(const double* const QPatch) const {
106 const idx4 idx_u(basisX, basisY, basisZ, numberOfVariables);
107 double Largest;
108 double pos[Dimensions];
109
110 for (int iz = 0; iz < basisZ; iz++)
111 for (int iy = 0; iy < basisY; iy++)
112 for (int ix = 0; ix < basisX; ix++) {
113 for(int id=0; id<Dimensions; id++) coordinate(pos, id);
114 double Lambda[numberOfVariables];
115
116 solver.eigenvalues(pos, QPatch + idx_u(iz,iy,ix,0), Lambda);
117 for(int k=0; k<numberOfVariables; k++) {
118 if(Lambda[k]>Largest) Largest = Lambda[k];
119 }
120 }
121
122 return Largest;
123 }
124
125 // Gets the largest value of the gaussLegendreNodes seperation. The attemp is slow and dumb.
126 double largest_dx() const {
127 double Largest;
128
129 for (int d = 0; d < Dimensions; d++)
130 for (int i = 0; i < basisSize; i++) {
131 double dx = sizeOfPatch[d] * SolverType::nodes[order][i];
132 if(dx > Largest) Largest = dx;
133 }
134
135 return Largest;
136 }
137
138 void Integrate(double time, double target_time, double* QPatch) const {
139 do {
140 // This timestep is hopefully very conservative
141 double dt = cfl_number / largest_eigenvalue(QPatch) / largest_dx();
142 dt = std::min(dt, target_time - time);
143
144 // RK3 for the time being:
145 double k[3][patchSize], QPatch1[patchSize], QPatch2[patchSize];
146
147 rhs(k[0], QPatch);
148 for(int i=0; i<patchSize; i++) QPatch1[i] = QPatch[i] + 0.5*dt*k[0][i];
149 rhs(k[1], QPatch1);
150 for(int i=0; i<patchSize; i++) QPatch2[i] = QPatch[i] - 1.0*dt*k[0][i] + 2.0*dt*k[1][i];
151 rhs(k[2], QPatch2);
152
153 for(int i=0; i<patchSize; i++) QPatch[i] += dt/6.0*(k[0][i] + 4.0*k[1][i] + 1.0*k[2][i]);
154
155 time += dt;
156 } while (time < target_time);
157 }
158}; // end of struct RungeKuttaIntegrator
159
160/* Other Runge Kutta schemes, could easily be supported:
161
162 CASE(1)
163 ! Explicit Euler
164 CALL Lh(k1,uh)
165 uh = uh + dt*k1
166 CASE(2)
167 ! RK2 scheme
168 CALL Lh(k1,uh)
169 CALL Lh(k2,uh+dt*k1)
170 uh = uh + dt/2.0*(k1+k2)
171 CASE(3)
172 ! Kutta's third order scheme
173 CALL Lh(k1,uh)
174 CALL Lh(k2,uh+0.5*dt*k1)
175 CALL Lh(k3,uh-1.0*dt*k1+2.0*dt*k2)
176 uh = uh + dt/6.0*( k1 + 4.0*k2 + 1.0*k3 )
177 CASE(4)
178 ! Classical RK4 scheme
179 CALL Lh(k1,uh)
180 CALL Lh(k2,uh+0.5*dt*k1)
181 CALL Lh(k3,uh+0.5*dt*k2)
182 CALL Lh(k4,uh+1.0*dt*k3)
183 uh = uh + dt/6.0*( k1 + 2.0*k2 + 2.0*k3 + k4 )
184*/
185
186
187} // c
188} // generic
189} // aderdg
190} // kernels
191
192#endif /* __KERNELS_STANDALONE_RK_INTEGRATOR_SVEN__ */
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx
-with-ExaHyPE2-benchmark
This file is part of the ExaHyPE project.
This is a small standalone RungeKutta Integrator suitable for integrating a single DG cell in time wi...
double largest_eigenvalue(const double *const QPatch) const
RungeKuttaIntegrator(const SolverType &solver, const dvec &offsetOfPatch, const dvec &sizeOfPatch, double cfl_number)
void coordinate(double *const pos, const int i) const
void Integrate(double time, double target_time, double *QPatch) const
void rhs(double *dudtPatch, const double *const QPatch) const