Loading...
Searching...
No Matches
16MSR_FOMSA.C
1#include "usmsrProblem.H"
2#include "ReducedUnsteadyMSR.H"
3#include "ITHACAstream.H"
4#include "LRSensitivity.H"
5#include "Tm.H"
6#include "Ptot.H"
7#include <Eigen/Dense>
8#include <math.h>
9#include <iomanip>
10#include <chrono>
11#include <ctime>
12#include "Tm_time.H"
13#include "Ptot_time.H"
14
15class msr : public usmsrProblem
16{
17 public:
18 explicit msr(int argc, char* argv[])
19 :
20 usmsrProblem(argc, argv),
21 U(_U()),
22 p(_p()),
23 flux(_flux()),
24 prec1(_prec1()),
25 prec2(_prec2()),
26 prec3(_prec3()),
27 prec4(_prec4()),
28 prec5(_prec5()),
29 prec6(_prec6()),
30 prec7(_prec7()),
31 prec8(_prec8()),
32 T(_T()),
33 dec1(_dec1()),
34 dec2(_dec2()),
35 dec3(_dec3()),
36 v(_v()),
37 D(_D()),
38 NSF(_NSF()),
39 A(_A()),
40 SP(_SP()),
41 TXS(_TXS())
42 {}
43
44 volVectorField& U;
45 volScalarField& p;
46 volScalarField& flux;
47 volScalarField& prec1;
48 volScalarField& prec2;
49 volScalarField& prec3;
50 volScalarField& prec4;
51 volScalarField& prec5;
52 volScalarField& prec6;
53 volScalarField& prec7;
54 volScalarField& prec8;
55 volScalarField& T;
56 volScalarField& dec1;
57 volScalarField& dec2;
58 volScalarField& dec3;
59 volScalarField& v;
60 volScalarField& D;
61 volScalarField& NSF;
62 volScalarField& A;
63 volScalarField& SP;
64 volScalarField& TXS;
65
66 ITHACAparameters* para = ITHACAparameters::getInstance(_mesh(), _runTime());
67 int NUproj = para->ITHACAdict->lookupOrDefault<int>("NUproj", 2);
68 int NPproj = para->ITHACAdict->lookupOrDefault<int>("NPproj", 2);
69 int NFluxproj = para->ITHACAdict->lookupOrDefault<int>("NFluxproj", 2);
70 int NPrecproj1 = para->ITHACAdict->lookupOrDefault<int>("NPrecproj1", 2);
71 int NPrecproj2 = para->ITHACAdict->lookupOrDefault<int>("NPrecproj2", 2);
72 int NPrecproj3 = para->ITHACAdict->lookupOrDefault<int>("NPrecproj3", 2);
73 int NPrecproj4 = para->ITHACAdict->lookupOrDefault<int>("NPrecproj4", 2);
74 int NPrecproj5 = para->ITHACAdict->lookupOrDefault<int>("NPrecproj5", 2);
75 int NPrecproj6 = para->ITHACAdict->lookupOrDefault<int>("NPrecproj6", 2);
76 int NPrecproj7 = para->ITHACAdict->lookupOrDefault<int>("NPrecproj7", 2);
77 int NPrecproj8 = para->ITHACAdict->lookupOrDefault<int>("NPrecproj8", 2);
78 int NTproj = para->ITHACAdict->lookupOrDefault<int>("NTproj", 2);
79 int NDecproj1 = para->ITHACAdict->lookupOrDefault<int>("NDecproj1", 2);
80 int NDecproj2 = para->ITHACAdict->lookupOrDefault<int>("NDecproj2", 2);
81 int NDecproj3 = para->ITHACAdict->lookupOrDefault<int>("NDecproj3", 2);
82 int NCproj = para->ITHACAdict->lookupOrDefault<int>("NCproj", 2);
83
84 double r1 = _beta1().value() / _betaTot().value();
85 double r2 = _beta2().value() / _betaTot().value();
86 double r3 = _beta3().value() / _betaTot().value();
87 double r4 = _beta4().value() / _betaTot().value();
88 double r5 = _beta5().value() / _betaTot().value();
89 double r6 = _beta6().value() / _betaTot().value();
90 double r7 = _beta7().value() / _betaTot().value();
91 double r8 = _beta8().value() / _betaTot().value();
92
93 void offlineSolve(std::string dir)
94 {
95 if (offline == false)
96 {
97 List<scalar> mu_now(3);
98 std::string folder = dir;
99
100 for (int i = 0; i < mu.cols(); i++)
101 {
102 folder.append(std::to_string(i));
103 folder.append("/");
104 _nu().value() = mu(0, i);
105 change_viscosity(mu(0, i));
106 _betaTot().value() = mu(1, i);
107 _beta1().value() = r1 * mu(1, i);
108 _beta2().value() = r2 * mu(1, i);
109 _beta3().value() = r3 * mu(1, i);
110 _beta4().value() = r4 * mu(1, i);
111 _beta5().value() = r5 * mu(1, i);
112 _beta6().value() = r6 * mu(1, i);
113 _beta7().value() = r7 * mu(1, i);
114 _beta8().value() = r8 * mu(1, i);
115 _decLam3().value() = mu(2, i);
116 mu_now[0] = mu(0, i);
117 mu_now[1] = mu(1, i);
118 mu_now[2] = mu(2, i);
119 truthSolve(mu_now, folder);
120 folder = dir;
121 restart();
122 }
123 }
124 else
125 {
126 readMSRfields(dir);
127 }
128 }
129
130};
131
132class Tmlocal : public Tm
133{
134 public:
135 explicit Tmlocal(int argc, char* argv[], int Nsampled)
136 :
137 Tm(argc, argv, Nsampled),
138 T(_T())
139 {}
140
141 volScalarField& T;
142};
143
144class Plocal: public Ptot
145{
146 public:
147 explicit Plocal(int argc, char* argv[], int Nsampled)
148 :
149 Ptot(argc, argv, Nsampled),
150 powerDens(_powerDens())
151 {}
152
153 volScalarField& powerDens;
154};
155
156int main(int argc, char* argv[])
157{
158 //object initialization
159 msr prova(argc, argv);
160 //inletIndex set
161 prova.inletIndex.resize(1, 2);
162 prova.inletIndex(0, 0) = 0;
163 prova.inletIndex(0, 1) = 0;
164 prova.Pnumber = 3;
165 prova.Tnumber = 1000;
166 prova.setParameters();
167 // define the distribution central and std. deviation values
168 double nu0 = 2.46E-06;
169 double betatot0 = 321.8E-05;
170 double dlam30 = 3.58e-04;
171 double signu = 0.1 / 3 * nu0;
172 double sigbeta = 0.1 / 3 * betatot0;
173 double sigdlam3 = 0.1 / 3 * dlam30;
174 // read range to adopt as training range
175 prova.mu_range = ITHACAstream::readMatrix("./range_mat.txt");
176 std::string dist = {"normal"};
177 // sample the values using ITHACA::sampling
178 prova.mu.row(0) = ITHACAsampling::samplingMC(dist, prova.mu_range(0, 0),
179 prova.mu_range(0, 1), nu0, signu, prova.Tnumber);
180 prova.mu.row(1) = ITHACAsampling::samplingMC(dist, prova.mu_range(1, 0),
181 prova.mu_range(1, 1), betatot0, sigbeta, prova.Tnumber);
182 prova.mu.row(2) = ITHACAsampling::samplingMC(dist, prova.mu_range(2, 0),
183 prova.mu_range(2, 1), dlam30, sigdlam3, prova.Tnumber);
184 double tstart = std::time(0);
185 Eigen::MatrixXd mu_f = prova.mu;
186
187 if (prova.offline == false)
188 {
189 ITHACAstream::exportMatrix(prova.mu, "mu_off", "eigen", "./ITHACAoutput");
190 }
191 else
192 {
193 mu_f = ITHACAstream::readMatrix("./ITHACAoutput/mu_off_mat.txt");
194 }
195
196 std::string folder = {"./ITHACAoutput/FOMoutput/"};
197 //perform the offline step
198 prova.offlineSolve(folder);
199 double tend = std::time(0);
200 double deltat_offline;
201 deltat_offline = difftime(tend, tstart);
202
203 if (prova.offline == false)
204 {
205 std::ofstream diff_t("./ITHACAoutput/t_SA=" + name(deltat_offline));
206 }
207
208 Info << "Offline stage time= " << deltat_offline << endl;
209 //sensitivity analysis
210 int Nparameters = prova.Pnumber;
211 int samplingPoints = prova.Tnumber;
212 LRSensitivity analisi(Nparameters, samplingPoints);
213 //load the sampling set
214 analisi.MatX.col(0) = mu_f.row(0);
215 analisi.MatX.col(1) = mu_f.row(1);
216 analisi.MatX.col(2) = mu_f.row(2);
217 analisi.getXstats();
218 //initialize Tm object
219 Tmlocal TmediaSA(argc, argv, samplingPoints);
220 Plocal PtotSA(argc, argv, samplingPoints);
221 //build the Model Output
222 TmediaSA.buildMO(folder);
223 PtotSA.buildMO(folder);
224 //assign it to VDSensitivity object
225 analisi.M = autoPtr<FofM>(& TmediaSA);
226 //load the model output in the VDSensitivity object
227 analisi.load_output();
228 analisi.getYstat();
229 Info << "Mean value and variance of the output:" << endl;
230 Info << analisi.Ey << "\t" << analisi.Vy << endl;
231 Info << "-------------" << endl;
232 analisi.getBetas();
233 Info << "analisi regression coefficients:" << endl;
234 Info << analisi.betas << endl;
235 Info << "-------------" << endl;
236 analisi.assessQuality();
237 Info << "quality: " << analisi.QI << endl;
238 Eigen::MatrixXd saveyout(samplingPoints, 1);
239 Eigen::MatrixXd saveymodel(samplingPoints, 1);
240 saveyout.col(0) = analisi.y;
241 saveymodel.col(0) = analisi.ylin;
242 ITHACAstream::exportMatrix(saveyout, "T_out", "eigen", "./ITHACAoutput");
243 ITHACAstream::exportMatrix(saveymodel, "T_model", "eigen", "./ITHACAoutput");
244 Eigen::MatrixXd coeffsT(4, 3);
245 coeffsT.setZero();
246 coeffsT(0, 0) = analisi.Ey;
247 coeffsT(0, 1) = analisi.Vy;
248 coeffsT(0, 2) = analisi.QI;
249 coeffsT.row(1) = analisi.EX;
250 coeffsT.row(2) = analisi.VX;
251 coeffsT.row(3) = analisi.betas;
252 ITHACAstream::exportMatrix(coeffsT, "coeffsT", "eigen", "./ITHACAoutput");
253 //assign it to VDSensitivity object
254 analisi.M = autoPtr<FofM>(& PtotSA);
255 //load the model output in the VDSensitivity object
256 analisi.load_output();
257 analisi.getYstat();
258 Info << "Mean value and variance of the output:" << endl;
259 Info << analisi.Ey << "\t" << analisi.Vy << endl;
260 Info << "-------------" << endl;
261 analisi.getBetas();
262 Info << "analisi regression coefficients:" << endl;
263 Info << analisi.betas << endl;
264 Info << "-------------" << endl;
265 analisi.assessQuality();
266 Info << "quality: " << analisi.QI << endl;
267 Eigen::MatrixXd saveyoutP(samplingPoints, 1);
268 Eigen::MatrixXd saveymodelP(samplingPoints, 1);
269 saveyoutP.col(0) = analisi.y;
270 saveymodelP.col(0) = analisi.ylin;
271 ITHACAstream::exportMatrix(saveyoutP, "P_out", "eigen", "./ITHACAoutput");
272 ITHACAstream::exportMatrix(saveymodelP, "P_model", "eigen", "./ITHACAoutput");
273 Eigen::MatrixXd coeffsP(4, 3);
274 coeffsP.setZero();
275 coeffsP(0, 0) = analisi.Ey;
276 coeffsP(0, 1) = analisi.Vy;
277 coeffsP(0, 2) = analisi.QI;
278 coeffsP.row(1) = analisi.EX;
279 coeffsP.row(2) = analisi.VX;
280 coeffsP.row(3) = analisi.betas;
281 ITHACAstream::exportMatrix(coeffsP, "coeffsP", "eigen", "./ITHACAoutput");
282 return 0;
283}
284
ITHACAstream.H
Header file of the ITHACAstream class, it contains the implementation of several methods for input ou...
ITHACAparameters::getInstance
static ITHACAparameters * getInstance()
Gets an instance of ITHACAparameters, to be used if the instance is already existing.
Definition ITHACAparameters.C:56
ITHACAsampling::samplingMC
static Eigen::VectorXd samplingMC(std::string pdftype, double &lowerE, double &upperE, double &distpara1, double &distpara2, label &Npoints)
Definition ITHACAsampling.C:5
Ptot
Definition Ptot.H:42
Ptot::_powerDens
autoPtr< volScalarField > _powerDens
List of pointers to power density field.
Definition Ptot.H:52
Tm
Definition Tm.H:42
Tm::_T
autoPtr< volScalarField > _T
List of pointers to temperature field.
Definition Tm.H:51
msrProblem::change_viscosity
void change_viscosity(double mu)
method to change the viscosity in UEqn
Definition msrProblem.C:1783
msrProblem::readMSRfields
void readMSRfields()
Method to read all the fieds of the MSR problem in the offline folder, it also reads mu_samples matri...
Definition msrProblem.C:1665
msrProblem::restart
void restart()
method to set all fields back to values in 0 folder
Definition msrProblem.C:2016
reductionProblem::offline
bool offline
Boolean variable, it is 1 if the Offline phase has already been computed, else 0.
Definition reductionProblem.H:99
reductionProblem::mu
Eigen::MatrixXd mu
Row matrix of parameters.
Definition reductionProblem.H:89
reductionProblem::truthSolve
void truthSolve()
Perform a TruthSolve.
Definition reductionProblem.C:97
usmsrProblem::usmsrProblem
usmsrProblem()
Construct Null.
Definition usmsrProblem.C:32
ITHACAstream::exportMatrix
void exportMatrix(Eigen::Matrix< T, -1, dim > &matrix, word Name, word type, word folder)
Export the reduced matrices in numpy (type=python), matlab (type=matlab) and txt (type=eigen) format ...
Definition ITHACAstream.C:55
ITHACAstream::readMatrix
List< Eigen::MatrixXd > readMatrix(word folder, word mat_name)
Read a three dimensional matrix from a txt file in Eigen format.
Definition ITHACAstream.C:372
Generated by   doxygen 1.16.1