NEDISS: Network Diffusion and Synchronization Simulator
GeneralSolver.h
1//
2// Created by m4zz31 on 1/11/21.
3//
4
5#ifndef CPPPROJCT_GENERALSOLVER_H
6#define CPPPROJCT_GENERALSOLVER_H
7
8#include <functional>
9
10
11#include "../GraphClasses/GeneralGraph.h"
12#include "../DifferentialEquations/GeneralDifferentialEquation.h"
13
14
15//
16// This file contains the definitions for the GeneralSolver Base Class
17// TODO: Document it, define Term1,2,3,4 so they can be overriden but are not mandatory for derived ones
18// the Document should include a mention that solvers do not support concurrent access :-) locks are not appropriate,
19// just copy the solver class in each thread. For GPU conversion, concurrent access can be allowed with some straightforward mods.
20
21typedef std::function<double(double, double, std::vector<double>, std::vector<double>, std::vector<double>)> ScalarFlow;
22typedef std::function<double(double a, std::vector<double> &b, std::vector<double> &c, std::vector<double> &d, Graph &g, ScalarFlow &F)> SolverOp;
23
24struct SolverConfig{
25 int s=0;
26 int d=0;
27 double P[4];
28 SolverConfig();
29};
30
31struct TimeStructure{
32 double h = 0.01;
33 double t = 0;
34};
35
36
37template <typename DIFFEQ, typename SOLVER>
38class GeneralSolver{
39public:
40 // Configuration specific to the solver structure :-)
41 double Params[4]; // up to 4 weights for runge-kutta
42 // for example, for Heun's method it would be {1/2,1/2,0,0}
43 // Note: in this example, k3 and k4 would still be computed if deg=4.
44 int deg; // the degree of the integration algorithm
45
46 std::string type;
47 bool requires_communication = false;
48 GeneralSolver();
49 GeneralSolver(std::string valtype);
50 GeneralSolver(std::string valtype, int d);
51 GeneralSolver(std::string valtype, int d, double * params);
52
53 // Instantiate the Differential Equation and prepare stuff if necessary
54 DIFFEQ DifferentialEquation;
55 SOLVER Solver;
56 void PostTemplateProcessing();
57
58
59 // main function: evolver
60 void evolve(double a,
61 std::vector<double> &b,
62 std::vector<double> &c,
63 std::vector<double> &d,
64 std::vector<double> &RK1,
65 std::vector<double> &RK2,
66 std::vector<double> &RK3,
67 std::vector<double> &RK4,
68 double &answer);
69 void Term1(double a,
70 std::vector<double> &b,
71 std::vector<double> &c,
72 std::vector<double> &d,
73 std::vector<double> &RK1);
74 void Term2(double a,
75 std::vector<double> &b,
76 std::vector<double> &c,
77 std::vector<double> &d,
78 std::vector<double> &RK1,
79 std::vector<double> &RK2);
80 void Term3(double a,
81 std::vector<double> &b,
82 std::vector<double> &c,
83 std::vector<double> &d,
84 std::vector<double> &RK1,
85 std::vector<double> &RK2,
86 std::vector<double> &RK3);
87 void Term4(double a,
88 std::vector<double> &b,
89 std::vector<double> &c,
90 std::vector<double> &d,
91 std::vector<double> &RK1,
92 std::vector<double> &RK2,
93 std::vector<double> &RK3,
94 std::vector<double> &RK4);
95
96
97 // Configuration specific to the time structure :-)
98 TimeStructure T;
99 void SetStep(double h);
100 void SetT0(double t0);
101 void EvolveTime();
102};
103
104// A flexible constructor that should be used carefully: d can be set up to the incorrect value
105template <typename DIFFEQ, typename SOLVER>
106GeneralSolver<DIFFEQ, SOLVER>::GeneralSolver(std::string valtype, int d, double * params){
107 type = valtype;
108 deg = d;
109 Params[0] = *(params+0);
110 Params[1] = *(params+1);
111 Params[2] = *(params+2);
112 Params[3] = *(params+3);
113 std::cout << "[WARNING] GeneralSolver constructor that explicitly uses four weights params assigns 'true' to 'requires_communication', change it if the current method does not require it. Current method: " << valtype << " with d=" << d << std::endl;
114 std::cout << std::flush;
115 requires_communication = DifferentialEquation.requiresCom(d);
116};
117
118// An easy constructor that yields Euler O(1) or RungeKutta O(2) ("Heun's Method")
119template <typename DIFFEQ, typename SOLVER>
120GeneralSolver<DIFFEQ, SOLVER>::GeneralSolver(std::string valtype){
121 // Types available should be:
122 // (1) Euler ('eu')
123 // (2) Runge-Kutta ('rk')
124 std::string methods_str = "(1) Euler ('eu'), (2) RungeKutta ('rk').";
125 type = std::move(valtype);
126 if (type == "rk") {
127 // Heun method's initialization :-)
128 deg = 2;
129 requires_communication = DifferentialEquation.requiresCom(deg);
130 Params[0] = 0.5;
131 Params[1] = 0.5;
132 Params[2] = 0;
133 Params[3] = 0;
134 } else if (type == "eu") {
135 deg = 1;
136 requires_communication = DifferentialEquation.requiresCom(deg);
137 Params[0] = 1;
138 Params[1] = 0;
139 Params[2] = 0;
140 Params[3] = 0;
141 } else {
142 error_report("Only support for Integration is available for the following methods:" + methods_str);
143 }
144};
145
146// The default constructor: Euler O(1)
147template <typename DIFFEQ, typename SOLVER>
148GeneralSolver<DIFFEQ, SOLVER>::GeneralSolver(){
149 type = "eu";
150 deg = 1;
151 Params[0] = 1;
152 Params[1] = 0;
153 Params[2] = 0;
154 Params[3] = 0;
155 requires_communication = DifferentialEquation.requiresCom(deg);
156}
157
158// "Easy" initializer that allows to receive only the name and the order.
159// For Euler it accepts any degree, for Runge Kutta only 'rk' and d=2 for Heun's Method.
160template <typename DIFFEQ, typename SOLVER>
161GeneralSolver<DIFFEQ, SOLVER>::GeneralSolver(std::string valtype, int d) {
162 // Types available should be:
163 // (1) Euler ('eu')
164 // (2) Runge-Kutta ('rk')
165 std::string methods_str = "(1) Euler ('eu'), (2) RungeKutta O(2) ('rk') a.k.a. Heun's Method";
166 type = std::move(valtype);
167 deg = d;
168 if (type == "rk") {
169 if (d == 2) {
170 // Heun method
171 Params[0] = 0.5;
172 Params[1] = 0.5;
173 Params[2] = 0;
174 Params[3] = 0;
175 requires_communication = DifferentialEquation.requiresCom(d);
176 } else {
177 printf("[FATAL] GeneralSolver(type='rk', int d) only accepts d=2 for Heun's method, but %d was the input.\n",d);
178 std::cout<<std::flush;
179 exit(1);
180 }
181 } else if (type == "eu") {
182 // Euler order 1 :-)
183 for (int i=0; i<4; ++i){
184 if (i+1 <= d){
185 Params[i] = 1;
186 } else {
187 Params[i] = 0;
188 }
189 requires_communication = DifferentialEquation.requiresCom(d);
190 }
191 if (d!=1) PRINTF_DBG("\n\n[WARNING]\n\nEuler od Order != 1 may not be available for some Equations, i.e. some equation classes may not have defined their field derivatives up to the requested order :O.\n\n");
192 } else {
193 error_report("Only support for Integration is available for the following methods:" + methods_str);
194 }
195};
196
197
198
199template <typename DIFFEQ, typename SOLVER>
200void GeneralSolver<DIFFEQ, SOLVER>::SetStep(double h){
201 T.h = h;
202}
203template <typename DIFFEQ, typename SOLVER>
204void GeneralSolver<DIFFEQ, SOLVER>::SetT0(double t0){
205 T.t = t0;
206}
207template <typename DIFFEQ, typename SOLVER>
208void GeneralSolver<DIFFEQ, SOLVER>::EvolveTime(){
209 T.t += T.h;
210}
211
212
213// The following are Class Functions that abstract away the call to the solver's
214// evolutions of the respective differential equations.
215template <typename DIFFEQ, typename SOLVER>
216void GeneralSolver<DIFFEQ, SOLVER>::evolve(double a,
217 std::vector<double> &b,
218 std::vector<double> &c,
219 std::vector<double> &d,
220 std::vector<double> &RK1,
221 std::vector<double> &RK2,
222 std::vector<double> &RK3,
223 std::vector<double> &RK4,
224 double &answer){
225 return Solver.evolve(T.t, T.h, a, b, c, d, DifferentialEquation,
226 DifferentialEquation.Specs, RK1, RK2, RK3, RK4, answer, &Params[0]);
227}
228
229
230template <typename DIFFEQ, typename SOLVER>
231void GeneralSolver<DIFFEQ, SOLVER>::Term1(double a,
232 std::vector<double> &b,
233 std::vector<double> &c,
234 std::vector<double> &d,
235 std::vector<double> &RK1){
236 return Solver.Term1(T.t, T.h, a, b, c, d, DifferentialEquation,
237 DifferentialEquation.Specs, RK1, &Params[0]);
238}
239
240
241template <typename DIFFEQ, typename SOLVER>
242void GeneralSolver<DIFFEQ, SOLVER>::Term2(double a,
243 std::vector<double> &b,
244 std::vector<double> &c,
245 std::vector<double> &d,
246 std::vector<double> &RK1,
247 std::vector<double> &RK2){
248 return Solver.Term2(T.t, T.h, a, b, c, d, DifferentialEquation,
249 DifferentialEquation.Specs, RK1, RK2, &Params[0]);
250}
251
252
253
254template <typename DIFFEQ, typename SOLVER>
255void GeneralSolver<DIFFEQ, SOLVER>::Term3(double a,
256 std::vector<double> &b,
257 std::vector<double> &c,
258 std::vector<double> &d,
259 std::vector<double> &RK1,
260 std::vector<double> &RK2,
261 std::vector<double> &RK3){
262 return Solver.Term3(T.t, T.h, a, b, c, d, DifferentialEquation,
263 DifferentialEquation.Specs, RK1, RK2, RK3, &Params[0]);
264}
265
266
267
268
269template <typename DIFFEQ, typename SOLVER>
270void GeneralSolver<DIFFEQ, SOLVER>::Term4(double a,
271 std::vector<double> &b,
272 std::vector<double> &c,
273 std::vector<double> &d,
274 std::vector<double> &RK1,
275 std::vector<double> &RK2,
276 std::vector<double> &RK3,
277 std::vector<double> &RK4){
278 return Solver.Term4(T.t, T.h, a, b, c, d, DifferentialEquation,
279 DifferentialEquation.Specs, RK1, RK2, RK3, RK4, &Params[0]);
280}
281
282
283
284template <typename DIFFEQ, typename SOLVER>
285void GeneralSolver<DIFFEQ, SOLVER>::PostTemplateProcessing(){
286 if (DifferentialEquation.RequiresBuilding){
287 DifferentialEquation.BuildForSolver();
288 };
289}
290
291
292
293
294#endif //CPPPROJCT_GENERALSOLVER_H
GeneralSolver
Definition: GeneralSolver.h:38
SolverConfig
Definition: GeneralSolver.h:24
TimeStructure
Definition: GeneralSolver.h:31
Generated by doxygen 1.9.2