1 #include <fstream>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <vector>
5 #include <string>
6 //MR CHANGE
7 #include <cstring>
8 //MR CHANGE END
9
10 #include "Input.h"
11 #include "common_c.h"
12
13 #include <FCMangle.h>
14 #define BC_setVars FortranCInterface_MODULE_(blowercontrol, bc_setvars, BLOWERCONTROL, BC_SETVARS)
15
16 using namespace std; //::cout;
17 void print_error_code(int ierr);
18 int SONFATH=0;
19 extern "C" char phasta_iotype[80];
20 //extern "C"
21 extern "C" void BC_setVars( int*, int*,
22 int*, int*,
23 double*, double*, double*,
24 double*, double*, double*,
25 double*, double*, double*,
26 double* );
27
input_fform(phSolver::Input & inp)28 int input_fform(phSolver::Input& inp)
29 {
30
31 int ierr = 0 ;
32 int i,j, n_tmp;
33
34 try {
35 if(workfc.myrank==workfc.master) {
36 printf("\n Complete Filename: %s \n", inp.GetDefaultFileName());
37 printf("\n Local Config: %s \n\n", inp.GetUserFileName());
38 }
39
40 #ifdef AMG
41
42 // AMG PARAMETERS
43
44 if ((string)inp.GetValue("Employ AMG") == "True" ) {
45
46 amgvari.irun_amg = 1;
47
48 amgvari.irun_amg_prec = inp.GetValue("Run AMG As CG-preconditioner");
49
50 amgvarr.strong_eps = inp.GetValue("Strong Criterion Eps");
51
52 amgvarr.ramg_eps = inp.GetValue("AMG Convergence Eps");
53
54 amgvarr.ramg_relax = inp.GetValue("AMG Relaxation Omega");
55 amgvarr.ramg_trunc = inp.GetValue("AMG Truncation Set");
56
57 amgvari.iamg_verb = inp.GetValue("AMG Verbosity");
58
59 amgvari.iamg_neg_sten = inp.GetValue("AMG Neg_Sten");
60
61 amgvari.iamg_nlevel = inp.GetValue("AMG Nlevel");
62
63 amgvari.iamg_c_solver = inp.GetValue("AMG Coarsest Solver");
64
65 amgvari.iamg_init = 0;
66 amgvari.iamg_setup_frez = inp.GetValue("AMG Freeze Setup");
67 if ((string)inp.GetValue("AMG Interpolation Type")=="Standard")
68 amgvari.iamg_interp = 1;
69 else
70 amgvari.iamg_interp = 0;
71 amgvari.maxnev = inp.GetValue("AMG GGB nev");
72 amgvari.maxncv = inp.GetValue("AMG GGB ncv");
73
74 if ((string)inp.GetValue("AMG Smoother Type")=="Gauss-Seidel")
75 amgvari.iamg_smoother = 1;
76 else if ((string)inp.GetValue("AMG Smoother Type")=="ChebyShev")
77 amgvari.iamg_smoother = 2;
78 else if ((string)inp.GetValue("AMG Smoother Type")=="MLS")
79 amgvari.iamg_smoother = 3;
80 amgvarr.ramg_chebyratio = inp.GetValue("AMG Chebyshev Eigenvalue ratio");
81
82 amgvari.mlsdeg = inp.GetValue("AMG MLS Degree");
83 amgvari.iamg_reduce = inp.GetValue("AMG Run Reduced Serial");
84 }
85 #endif
86
87 /////////////////////////////chen Sep 25 2009 Flow Control Parameters ////////
88 // Take BC from IC at inlet
89 ctrlvari.iI2Binlet = (int)inp.GetValue("Take BC from IC at Inlet");
90 if(ctrlvari.iI2Binlet > 0 )
91 {
92 ctrlvar.inletVelX = inp.GetValue("Inlet Bulk x Velocity");
93 }
94 // set uniform outlet pressure
95 // ctrlvari.isetOutPres = (int)inp.GetValue("Set Outlet Pressure");
96 // if(ctrlvari.isetOutPres > 0 )
97
98 ductvari.isetOutletID = (int)inp.GetValue("Duct Outlet ID");
99 if(ductvari.isetOutletID > 0 )
100 ctrlvar.outPres1 = inp.GetValue("Duct Uniform Outlet Pressure");
101
102 // Override Eddy Vicosity IC and BC
103 ductvari.isetEV_IC_BC=(int)inp.GetValue("Override Eddy Viscosity");
104 if(ductvari.isetEV_IC_BC==1){
105 ductvar.evis_IC_BC=inp.GetValue("Eddy Viscosity Value for Override");
106 }
107
108 if(ductvari.isetEVramp = (int)inp.GetValue("Specify Initial Eddy Viscosity Ramp")){
109 ductvar.EVrampXmin =inp.GetValue("Initial Scalar 1 ramp start");
110 ductvar.EVrampXmax =inp.GetValue("Initial Scalar 1 ramp end");
111 ductvar.EVrampMax =inp.GetValue("Initial Scalar 1 high");
112 ductvar.EVrampMin =inp.GetValue("Initial Scalar 1 low");
113 }
114
115 // set initial condition
116 ctrlvari.isetInitial=(int)inp.GetValue("Specify Initial Conditions");
117 if(ctrlvari.isetInitial==1){
118 ctrlvar.xvel_ini = inp.GetValue("Initial X Velocity");
119 ctrlvar.yvel_ini = inp.GetValue("Initial Y Velocity");
120 ctrlvar.zvel_ini = inp.GetValue("Initial Z Velocity");
121 ctrlvar.temp_ini = inp.GetValue("Initial Temp");
122 ctrlvar.pres_ini = inp.GetValue("Initial Pressure");
123 ctrlvar.evis_ini = inp.GetValue("Initial Scalar 1");
124 }
125
126
127
128
129 //initial condition for duct
130 ductvari.isetInitial_Duct=(int)inp.GetValue("Set Initial Condition for Duct");
131 //inlet condition for duct
132 ductvari.isetInlet_Duct=(int)inp.GetValue("Set Inlet Condition for Duct");
133
134 //surfID, t_cycle, t_riseTime, t_fallTime, t_fullOn, vmax, vmin, T, nu, deltaBL, enable
135 n_tmp = (int) inp.GetValue("Number of Blower Surfaces"); //BC_setNBlower(&n_tmp);
136
137 if(n_tmp > 0){
138 vector<int> *ivec[3];
139 vector<double> *dvec[10];
140
141 for(i = 0; i < 3; i++) ivec[i] = new vector<int>;
142 for(i = 0; i < 10; i++) dvec[i] = new vector<double>;
143
144 *ivec[0] = inp.GetValue("Blower Mode");
145 *ivec[1] = inp.GetValue("Blower Surface ID");
146 *ivec[2] = inp.GetValue("Blower Enable");
147
148 *dvec[0] = inp.GetValue("Blower Cycle Period");
149 *dvec[1] = inp.GetValue("Blower Full On Period");
150 *dvec[2] = inp.GetValue("Blower Rise Time");
151 *dvec[3] = inp.GetValue("Blower Fall Time");
152 *dvec[4] = inp.GetValue("Blower Maximum u_normal");
153 *dvec[5] = inp.GetValue("Blower Minimum u_normal");
154 *dvec[6] = inp.GetValue("Blower Temperature");
155 *dvec[7] = inp.GetValue("Blower Eddy Viscosity");
156 *dvec[8] = inp.GetValue("Blower BL Thickness");
157 *dvec[9] = inp.GetValue("Blower BL Thickness (scalar)");
158
159 BC_setVars( &n_tmp,
160 &(*ivec[0])[0], //mode
161 &(*ivec[1])[0], //surfID
162 &(*ivec[2])[0], //enable
163 &(*dvec[0])[0], //t_cycle
164 &(*dvec[1])[0], //t_fullOn
165 &(*dvec[2])[0], //t_riseTime
166 &(*dvec[3])[0], //t_fallTime
167 &(*dvec[4])[0], //vmax
168 &(*dvec[5])[0], //vmin
169 &(*dvec[6])[0], //T
170 &(*dvec[7])[0], //nu
171 &(*dvec[8])[0], //delta BL velocity
172 &(*dvec[9])[0] ); //delta BL scalar
173
174 }
175
176 //suction condition for duct
177 ductvari.isetSuctionID_Duct=(int)inp.GetValue("Duct Set Suction Surface ID"); //suction patch surface IDs
178 if(ductvari.isetSuctionID_Duct > 0){
179 ductvari.suctionVbottom = inp.GetValue("Duct Bottom Suction Normal Velocity");
180 ductvari.suctionVside_lower = inp.GetValue("Duct Lower Side Suction Normal Velocity");
181 ductvari.suctionVside_upper = inp.GetValue("Duct Upper Side Surface Normal Velocity");
182 ductvari.suctionVtop = inp.GetValue("Duct Top Surface Normal Velocity");
183 }
184
185 // duct geometry type
186 ductvari.iDuctgeometryType = (int)inp.GetValue("Duct Geometry Type");
187
188 ///////////////////////////////////////////////////////////////////////////////
189
190
191 // Disabled Features
192
193 conpar.iALE = inp.GetValue("iALE");
194 conpar.icoord = inp.GetValue("icoord");
195 conpar.irs = inp.GetValue("irs");
196 conpar.iexec = inp.GetValue("iexec");
197 timpar.ntseq = inp.GetValue("ntseq");
198 solpar.imap = inp.GetValue("imap");
199
200
201 // Solution Control Keywords
202
203 if((string)inp.GetValue("Equation of State") == "Incompressible") matdat.matflg[0][0] =-1 ;
204 if((string)inp.GetValue("Equation of State") == "Compressible") matdat.matflg[0][0] =0;
205 inpdat.Delt[0] = inp.GetValue("Time Step Size");
206 inpdat.nstep[0] = inp.GetValue("Number of Timesteps");
207 if((string)inp.GetValue("Viscous Control")=="Viscous") conpar.navier=1 ; else conpar.navier=0;
208
209 if ((string)inp.GetValue("Turbulence Model") == "No-Model" ) {
210 turbvari.irans = 0;
211 turbvari.iles = 0;
212 } else if ((string)inp.GetValue("Turbulence Model") == "LES" ) {
213 turbvari.iles = 1;
214 turbvari.irans = 0;
215 } else if ((string)inp.GetValue("Turbulence Model") == "RANS-SA" ) {
216 turbvari.iles = 0;
217 turbvari.irans = -1;
218 } else if ((string)inp.GetValue("Turbulence Model") == "RANS" ) {
219 turbvari.iles = 0;
220 turbvari.irans = -1; // assume S-A for backward compatibility
221 } else if ((string)inp.GetValue("Turbulence Model") == "RANS-KE" ) {
222 turbvari.iles = 0;
223 turbvari.irans = -2;
224 } else if ((string)inp.GetValue("Turbulence Model") == "DES97" ) {
225 turbvari.iles = -1;
226 turbvari.irans = -1;
227 } else if ((string)inp.GetValue("Turbulence Model") == "DDES" ) {
228 turbvari.iles = -2;
229 turbvari.irans = -1;
230
231 } else {
232 cout << " Turbulence Model: Only Legal Values ( No-Model, LES, RANS-SA, RANS-KE, DES97, DDES )";
233 cout << endl;
234 exit(1);
235 }
236
237 // if (turbvari.iles*turbvari.irans!=0) turbvar.eles=
238 // inp.GetValue("DES Edge Length");
239
240 if (turbvari.irans<0 && turbvari.iles<0)
241 turbvar.DES_SA_hmin=(double)inp.GetValue("DES SA Minimum Edge Length");
242
243 int solflow, solheat , solscalr, ilset;
244 ((string)inp.GetValue("Solve Flow") == "True")? solflow=1:solflow=0;
245 ((string)inp.GetValue("Solve Heat") == "True")? solheat=1:solheat=0;
246 //for compressible solheat= False so
247 if((string)inp.GetValue("Equation of State") == "Compressible") solheat=0;
248 ilset = (int)inp.GetValue("Solve Level Set");
249 solscalr = (int)inp.GetValue("Solve Scalars");
250 solscalr += ilset;
251 if(turbvari.irans == -1) solscalr++;
252 if(turbvari.irans == -2) solscalr=solscalr+2;
253 if ( solscalr > 4 ) {
254 cout << " Only Four Scalars are supported \n";
255 cout <<" Please reduce number of scalars \n";
256 exit(1);
257 }
258 inpdat.impl[0] = 10*solflow+solscalr*100+solheat;
259
260 levlset.iLSet = ilset;
261 if( ilset > 0) {
262 levlset.epsilon_ls = inp.GetValue("Number of Elements Across Interface");
263 levlset.epsilon_lsd = inp.GetValue("Number of Elements Across Interface for Redistancing");
264 levlset.dtlset = inp.GetValue("Pseudo Time step for Redistancing");
265 levlset.iExpLSSclr2 = inp.GetValue("Explicit Solve for Redistance Field");
266 levlset.iExpLSSclr1 = inp.GetValue("Explicit Solve for Scalar 1 Field");
267 if ((string)inp.GetValue("Apply Volume Constraint") == "True" ) {
268 levlset.ivconstraint = 1; }
269 else if((string)inp.GetValue("Apply Volume Constraint") == "False" ) {
270 levlset.ivconstraint = 0; }
271 else {
272 cout << "Apply Volume Constraint: Only Legal Values (True, False) ";
273 cout << endl;
274 exit(1);
275 }
276 }
277
278 vector<double> vec;
279
280 // OUTPUT CONTROL KEY WORDS.
281
282 conpar.necho = inp.GetValue("Verbosity Level");
283 outpar.ntout = inp.GetValue("Number of Timesteps between Restarts");
284 outpar.nsynciofiles = inp.GetValue("Number of SyncIO Files");
285 if((string)inp.GetValue("Print Statistics") == "True") outpar.ioform = 2;
286 else outpar.ioform = 1;
287
288 if((string)inp.GetValue("Print Wall Fluxes") == "True") outpar.iowflux = 1;
289 else outpar.iowflux = 0;
290
291 if((string)inp.GetValue("Print FieldView") == "True") outpar.iofieldv = 1;
292 else outpar.iofieldv = 0;
293
294 if((string)inp.GetValue("Print ybar") == "True") outpar.ioybar = 1;
295 else outpar.ioybar = 0;
296
297 if((string)inp.GetValue("Print vorticity") == "True") outpar.ivort = 1;
298 else outpar.ivort = 0;
299
300 outpar.nstepsincycle = inp.GetValue("Number of Steps in a Cycle");
301 outpar.nphasesincycle = inp.GetValue("Number of Phases in a Cycle");
302 outpar.ncycles_startphaseavg = inp.GetValue("Number of Initial Cycles to Skip in Phase Average");
303
304 strcpy( outpar.iotype , ((string)inp.GetValue("Data Block Format")).c_str());
305 strcpy( phasta_iotype , ((string)inp.GetValue("Data Block Format")).c_str());
306 SONFATH = inp.GetValue("Number of Father Nodes");
307
308 if((string)inp.GetValue("Print Residual at End of Step") == "True") genpar.lstres = 1;
309 else genpar.lstres = 0;
310
311 if((string)inp.GetValue("Print Error Indicators") == "True") turbvar.ierrcalc = 1;
312 else turbvar.ierrcalc = 0;
313
314 if((string)inp.GetValue("Print Velocity Hessian") == "True") turbvar.ihessian = 1;
315 else turbvar.ihessian = 0;
316
317 if ( turbvar.ierrcalc == 1 )
318 turbvari.ierrsmooth = inp.GetValue("Number of Error Smoothing Iterations");
319
320 for(i=0;i<MAXSURF+1; i++) aerfrc.nsrflist[i] = 0;
321 int nsrfCM = inp.GetValue("Number of Force Surfaces");
322 if (nsrfCM > 0) {
323 vector<int> ivec = inp.GetValue("Surface ID's for Force Calculation");
324 for(i=0; i< nsrfCM; i++){
325 aerfrc.nsrflist[ivec[i]] = 1;
326 // cout <<"surface in force list "<< ivec[i] << endl;
327 }
328 ivec.erase(ivec.begin(),ivec.end());
329 }
330
331 aerfrc.isrfIM = inp.GetValue("Surface ID for Integrated Mass");
332
333 outpar.iv_rankpercore = inp.GetValue("Ranks per core");
334 outpar.iv_corepernode = inp.GetValue("Cores per node");
335
336 turbvari.iramp=0;
337 if((string)inp.GetValue("Ramp Inflow") == "True") turbvari.iramp=1;
338 if(turbvari.iramp == 1) {
339 vec = inp.GetValue("Mdot Ramp Inflow Start and Stop");
340 for(i=0; i<2 ; i++){
341 turbvar.rampmdot[0][i]=vec[i];
342 }
343 vec = inp.GetValue("Mdot Ramp Lower FC Start and Stop");
344 for(i=0; i<2 ; i++){
345 turbvar.rampmdot[1][i]=vec[i];
346 }
347 vec = inp.GetValue("Mdot Ramp Upper FC Start and Stop");
348 for(i=0; i<2 ; i++){
349 turbvar.rampmdot[2][i]=vec[i];
350 }
351 }
352
353 //Limiting
354 vec = inp.GetValue("Limit u1");
355 for(i=0; i<3 ; i++){
356 turbvar.ylimit[0][i] = vec[i];
357 }
358 vec.erase(vec.begin(),vec.end());
359
360 vec = inp.GetValue("Limit u2");
361 for(i=0; i<3 ; i++){
362 turbvar.ylimit[1][i] = vec[i];
363 }
364 vec.erase(vec.begin(),vec.end());
365
366 vec = inp.GetValue("Limit u3");
367 for(i=0; i<3 ; i++){
368 turbvar.ylimit[2][i] = vec[i];
369 }
370 vec.erase(vec.begin(),vec.end());
371
372 vec = inp.GetValue("Limit Pressure");
373 for(i=0; i<3 ; i++){
374 turbvar.ylimit[3][i] = vec[i];
375 }
376 vec.erase(vec.begin(),vec.end());
377
378 vec = inp.GetValue("Limit Temperature");
379 for(i=0; i<3 ; i++){
380 turbvar.ylimit[4][i] = vec[i];
381 }
382 vec.erase(vec.begin(),vec.end());
383
384 //Material Properties Keywords
385 matdat.nummat = levlset.iLSet+1;
386 if((string)inp.GetValue("Shear Law") == "Constant Viscosity")
387 for(i=0; i < levlset.iLSet+1; i++) matdat.matflg[i][1] = 0;
388
389 if((string)inp.GetValue("Bulk Viscosity Law") == "Constant Bulk Viscosity")
390 for(i=0; i < levlset.iLSet+1; i++) matdat.matflg[i][2] = 0;
391
392 mmatpar.pr = inp.GetValue("Prandtl Number");
393
394 if((string)inp.GetValue("Conductivity Law") == "Constant Conductivity")
395 for(i=0; i < levlset.iLSet+1; i++) matdat.matflg[i][3] = 0;
396
397 vec = inp.GetValue("Density");
398 for(i=0; i< levlset.iLSet +1 ; i++){
399 matdat.datmat[i][0][0] = vec[i];
400 }
401 vec.erase(vec.begin(),vec.end());
402
403 vec = inp.GetValue("Viscosity");
404 for(i=0; i< levlset.iLSet +1 ; i++){
405 matdat.datmat[i][1][0] = vec[i];
406 }
407 vec.erase(vec.begin(),vec.end());
408
409 // vec = inp.GetValue("Specific Heat");
410 for(i=0; i< levlset.iLSet +1 ; i++){
411 matdat.datmat[i][2][0] = 0;
412 }
413 // vec.erase(vec.begin(),vec.end());
414
415 vec = inp.GetValue("Thermal Conductivity");
416 for(i=0; i< levlset.iLSet +1 ; i++){
417 matdat.datmat[i][3][0] = vec[i];
418 }
419 vec.erase(vec.begin(),vec.end());
420
421 vec = inp.GetValue("Scalar Diffusivity");
422 for(i=0; i< solscalr ; i++){
423 sclrs.scdiff[i] = vec[i];
424 }
425 vec.erase(vec.begin(),vec.end());
426
427 if((string)inp.GetValue("Zero Mean Pressure") == "True")
428 turbvar.pzero=1;
429
430 turbvar.rmutarget = inp.GetValue("Target Viscosity For Step NSTEP");
431
432 if ( (string)inp.GetValue("Body Force Option") == "None" ) {
433 for( i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 0;
434 }
435 else if ( (string)inp.GetValue("Body Force Option") == "Vector" ) {
436 for( i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 1;
437 }
438 else if ( (string)inp.GetValue("Body Force Option") == "User e3source.f" ) {
439 for( i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 3;
440 }
441 else if ( (string)inp.GetValue("Body Force Option") == "Boussinesq" ) {
442 for(i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 2;
443 }
444 else if ( (string)inp.GetValue("Body Force Option") == "Cooling Analytic" ) {
445 for(i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 4;
446 ctrlvar.xvel_ini = inp.GetValue("Initial X Velocity");
447 ctrlvar.yvel_ini = inp.GetValue("Initial Y Velocity");
448 ctrlvar.zvel_ini = inp.GetValue("Initial Z Velocity");
449 ctrlvar.temp_ini = inp.GetValue("Initial Temp");
450 ctrlvar.pres_ini = inp.GetValue("Initial Pressure");
451 }
452 else if ( (string)inp.GetValue("Body Force Option") == "Cooling Initial Condition" ) {
453 for(i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 5;
454 }
455
456 // the following block of stuff is common to all cooling type sponges.
457 // Specific stuff belongs in the conditionals above
458
459 if(matdat.matflg[0][4] >=4) {
460 spongevar.betamax = inp.GetValue("Maximum Value of Sponge Parameter");
461 spongevar.zinsponge = inp.GetValue("Inflow Cooling Sponge Ends at z");
462 spongevar.zoutsponge= inp.GetValue("Outflow Cooling Sponge Begins at z");
463 spongevar.radsponge = inp.GetValue("Radial Cooling Sponge Begins at r");
464 spongevar.grthosponge = inp.GetValue("Sponge Growth Coefficient Outflow");
465 spongevar.grthisponge = inp.GetValue("Sponge Growth Coefficient Inflow");
466
467
468 spongevar.spongecontinuity = 0;
469 spongevar.spongemomentum1 = 0;
470 spongevar.spongemomentum2 = 0;
471 spongevar.spongemomentum3 = 0;
472 spongevar.spongeenergy = 0;
473
474 if((string)inp.GetValue("Sponge for Continuity Equation") == "True")
475 spongevar.spongecontinuity = 1;
476 if((string)inp.GetValue("Sponge for x Momentum Equation") == "True")
477 spongevar.spongemomentum1 = 1;
478 if((string)inp.GetValue("Sponge for y Momentum Equation") == "True")
479 spongevar.spongemomentum2 = 1;
480 if((string)inp.GetValue("Sponge for z Momentum Equation") == "True")
481 spongevar.spongemomentum3 = 1;
482 if((string)inp.GetValue("Sponge for Energy Equation") == "True")
483 spongevar.spongeenergy = 1;
484
485 }
486
487 vec = inp.GetValue("Body Force");
488 for(i=0; i< levlset.iLSet +1 ; i++){
489 matdat.datmat[i][4][0] = vec[0+i*3];
490 matdat.datmat[i][4][1] = vec[1+i*3];
491 matdat.datmat[i][4][2] = vec[2+i*3];
492 }
493 vec.erase(vec.begin(),vec.end());
494
495 vec = inp.GetValue("Body Force Pressure Gradient");
496 for(i=0; i< levlset.iLSet +1 ; i++){
497 matdat.datmat[i][6][0] = vec[0+i*3];
498 matdat.datmat[i][6][1] = vec[1+i*3];
499 matdat.datmat[i][6][2] = vec[2+i*3];
500 }
501 vec.erase(vec.begin(),vec.end());
502
503 if ( (string)inp.GetValue("Surface Tension Option") == "No" ){
504 genpar.isurf = 0;
505 }
506 else if ((string)inp.GetValue("Surface Tension Option") == "Yes" ){
507 genpar.isurf = 1;
508 }
509 else {
510 cout << " Surface Tension: Only Legal Values (Yes, No) ";
511 cout << endl;
512 exit(1);
513 }
514 if( genpar.isurf > 0) {
515 genpar.Bo = inp.GetValue("Bond Number");
516 }
517
518 genpar.EntropyPressure = inp.GetValue("Entropy Form of Pressure Constraint on Weight Space");
519
520
521 if ( (string)inp.GetValue("Rotating Frame of Reference") == "True" ) {
522 matdat.matflg[0][5] = 1;
523 vec = inp.GetValue("Rotating Frame of Reference Rotation Rate");
524 matdat.datmat[0][5][0] = vec[0];
525 matdat.datmat[0][5][1] = vec[1];
526 matdat.datmat[0][5][2] = vec[2];
527 vec.erase(vec.begin(),vec.end());
528 }
529 else {
530 matdat.matflg[0][5] = 0;
531 matdat.datmat[0][5][0] = 0.;
532 matdat.datmat[0][5][1] = 0.;
533 matdat.datmat[0][5][2] = 0.;
534 }
535
536
537 //Linear Solver parameters
538 conpar.usingpetsc=0; // default is to have PETSc off
539 incomp.iprjFlag = 0; incomp.ipresPrjFlag=0; inpdat.svLSFlag=0;
540 inpdat.leslib=0;
541 if( (string)inp.GetValue("Solver Type") =="svLS" ){
542 inpdat.svLSFlag = 1; }
543 if( (string)inp.GetValue("Solver Type") =="ACUSIM with P Projection" ){
544 inpdat.leslib=1; incomp.iprjFlag = 0; incomp.ipresPrjFlag=1;}
545 else if ( (string)inp.GetValue("Solver Type") =="ACUSIM" ){
546 inpdat.leslib=1; incomp.iprjFlag = 0; incomp.ipresPrjFlag=0;}
547 else if( (string)inp.GetValue("Solver Type") =="ACUSIM with Velocity Projection" ){
548 inpdat.leslib=1; incomp.iprjFlag = 1; incomp.ipresPrjFlag=0;}
549 else if( (string)inp.GetValue("Solver Type") =="ACUSIM with Full Projection" ){
550 inpdat.leslib=1; incomp.iprjFlag = 1; incomp.ipresPrjFlag=1;}
551 else if( (string)inp.GetValue("Solver Type") =="GMRES Matrix Free"){
552 inpdat.impl[0] += 10*solflow;}
553 else if( (string)inp.GetValue("Solver Type") =="GMRES EBE"){
554 inpdat.impl[0] += 20*solflow;}
555 else if( (string)inp.GetValue("Solver Type") =="PETSc"){
556 conpar.usingpetsc=1;}
557 //GMRES sparse is assumed default and has the value of 10, MFG 20,
558 // EBE 30
559
560
561 // inpdat.niter[0] = inp.GetValue("Number of Solves per Time Step");
562 solpar.nGMRES = inp.GetValue("Number of GMRES Sweeps per Solve");
563 solpar.Kspace = inp.GetValue("Number of Krylov Vectors per GMRES Sweep");
564 inpdat.LHSupd[0] = inp.GetValue("Number of Solves per Left-hand-side Formation");
565 inpdat.epstol[0] = inp.GetValue("Tolerance on Momentum Equations");
566 incomp.prestol = inp.GetValue("Tolerance on ACUSIM Pressure Projection");
567 incomp.minIters = inp.GetValue("Minimum Number of Iterations per Nonlinear Iteration");
568 incomp.maxIters = inp.GetValue("Maximum Number of Iterations per Nonlinear Iteration");
569 inpdat.deltol[0][0]=inp.GetValue("Velocity Delta Ratio");
570 inpdat.deltol[1][0]=inp.GetValue("Pressure Delta Ratio");
571 incomp.nPrjs = inp.GetValue("Number of Velocity Projection Vectors");
572 incomp.nPresPrjs = inp.GetValue("Number of Pressure Projection Vectors");
573 incomp.iverbose = inp.GetValue("ACUSIM Verbosity Level");
574
575 if(solheat==1){
576 inpdat.epstol[1]=inp.GetValue("Temperature Solver Tolerance");
577 inpdat.LHSupd[1]=inp.GetValue("Number of Solves of Temperature per Left-hand-side Formation");
578 }
579
580 // The following is where you should put any inputs that are able to
581 // input differently for each scalar. It is a little tedious in the code
582 // but it should make the solver.inp easier to understand. Note this will
583 // require some care with regression tests.
584
585
586 if(solscalr>0){
587 inpdat.epstol[2]=inp.GetValue("Scalar 1 Solver Tolerance");
588 inpdat.LHSupd[2]=inp.GetValue("Number of Solves of Scalar 1 per Left-hand-side Formation");
589
590 vec = inp.GetValue("Limit Scalar 1");
591 for(i=0; i<3 ; i++){
592 turbvar.ylimit[5][i] = vec[i];
593 }
594 vec.erase(vec.begin(),vec.end());
595 }
596
597 if(solscalr>1){
598 inpdat.epstol[3]=inp.GetValue("Scalar 2 Solver Tolerance");
599 inpdat.LHSupd[3]=inp.GetValue("Number of Solves of Scalar 2 per Left-hand-side Formation");
600
601 vec = inp.GetValue("Limit Scalar 2");
602 for(i=0; i<3 ; i++){
603 turbvar.ylimit[6][i] = vec[i];
604 }
605 vec.erase(vec.begin(),vec.end());
606 }
607
608 if(solscalr>2){
609 inpdat.epstol[4]=inp.GetValue("Scalar 3 Solver Tolerance");
610 inpdat.LHSupd[4]=inp.GetValue("Number of Solves of Scalar 3 per Left-hand-side Formation");
611
612 vec = inp.GetValue("Limit Scalar 3");
613 for(i=0; i<3 ; i++){
614 turbvar.ylimit[7][i] = vec[i];
615 }
616 vec.erase(vec.begin(),vec.end());
617 }
618
619 if(solscalr>3){
620 inpdat.epstol[5]=inp.GetValue("Scalar 4 Solver Tolerance");
621 inpdat.LHSupd[5]=inp.GetValue("Number of Solves of Scalar 4 per Left-hand-side Formation");
622
623 vec = inp.GetValue("Limit Scalar 4");
624 for(i=0; i<3 ; i++){
625 turbvar.ylimit[8][i] = vec[i];
626 }
627 vec.erase(vec.begin(),vec.end());
628 }
629
630 // DISCRETIZATION CONTROL
631
632 genpar.ipord = inp.GetValue("Basis Function Order");
633 if((string)inp.GetValue("Time Integration Rule") == "First Order")
634 inpdat.rhoinf[0] = -1 ;
635 else inpdat.rhoinf[0] = (double)inp.GetValue("Time Integration Rho Infinity");
636 if((string)inp.GetValue("Predictor at Start of Step")=="Same Velocity")
637 genpar.ipred = 1;
638 if((string)inp.GetValue("Predictor at Start of Step")=="Zero Acceleration")
639 genpar.ipred = 2;
640 if((string)inp.GetValue("Predictor at Start of Step")=="Same Acceleration")
641 genpar.ipred = 3;
642 if((string)inp.GetValue("Predictor at Start of Step")=="Same Delta")
643 genpar.ipred = 4;
644
645 if((string)inp.GetValue("Weak Form") == "Galerkin")
646 solpar.ivart = 1;
647 if((string)inp.GetValue("Weak Form") == "SUPG")
648 solpar.ivart = 2;
649
650 if((string)inp.GetValue("Flow Advection Form") == "Convective")
651 solpar.iconvflow = 2;
652 else if((string)inp.GetValue("Flow Advection Form") == "Conservative")
653 solpar.iconvflow = 1;
654 if((string)inp.GetValue("Scalar Advection Form") == "Convective")
655 solpar.iconvsclr = 2;
656 else if((string)inp.GetValue("Scalar Advection Form") == "Conservative")
657 solpar.iconvsclr = 1;
658 if((string)inp.GetValue("Use Conservative Scalar Convection Velocity") == "True")
659 sclrs.consrv_sclr_conv_vel = 1;
660 else if((string)inp.GetValue("Use Conservative Scalar Convection Velocity") == "False")
661 sclrs.consrv_sclr_conv_vel = 0;
662 // TAU INPUT
663 if((string)inp.GetValue("Tau Matrix") == "Diagonal-Shakib")
664 genpar.itau = 0;
665 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Franca")
666 genpar.itau =1;
667 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Jansen(dev)")
668 genpar.itau = 2;
669 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Compressible")
670 genpar.itau = 3;
671 else if((string)inp.GetValue("Tau Matrix") == "Matrix-Mallet")
672 genpar.itau = 10;
673 else if((string)inp.GetValue("Tau Matrix") == "Matrix-Modal")
674 genpar.itau = 11;
675
676 genpar.dtsfct = inp.GetValue("Tau Time Constant");
677 genpar.taucfct = inp.GetValue("Tau C Scale Factor");
678
679 genpar.iLHScond = inp.GetValue("LHS BC heat flux enable");
680
681 // FLOW DISCONTINUITY CAPTURING
682
683 if((string)inp.GetValue("Discontinuity Capturing") == "Off") solpar.iDC = 0;
684 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-mallet") solpar.iDC = 1;
685 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-quadratic") solpar.iDC = 2;
686 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-minimum") solpar.iDC = 3;
687 else {
688 cout<< "Condition not defined for Discontinuity Capturing \n ";
689 exit(1);
690 }
691
692 // SCALAR DISCONTINUITY CAPTURING
693
694 vector<int> ivec = inp.GetValue("Scalar Discontinuity Capturing");
695 for(i=0; i< 2; i++) solpar.idcsclr[i] = ivec[i];
696 ivec.erase(ivec.begin(),ivec.end());
697
698
699 // if((string)inp.GetValue("Scalar Discontinuity Capturing") == "No") solpar.idcsclr = 0;
700 // else if((string)inp.GetValue("Scalar Discontinuity Capturing") == "1") solpar.idcsclr = 1;
701 // else if((string)inp.GetValue("Scalar Discontinuity Capturing") == "2") solpar.idcsclr = 2;
702 // else {
703 // cout<< "Condition not defined for Scalar Discontinuity Capturing \n ";
704 // exit(1);
705 // }
706 if((string)inp.GetValue("Include Viscous Correction in Stabilization") == "True")
707 {
708 if(genpar.ipord == 1 ) genpar.idiff = 1;
709 else genpar.idiff = 2;
710 }
711 else { genpar.idiff = 0;}
712
713 // ------ Only For duct S duct Project ------------------------------------------------
714 genpar.irampViscOutlet = (int)inp.GetValue("Ramp Up Viscosity Near Outlet");
715
716 genpar.istretchOutlet = (int)inp.GetValue("Stretch X Coordinate Near Outlet");
717 // -----------------------------------------------------------------------------------
718
719 genpar.iremoveStabTimeTerm = (int)inp.GetValue("Remove Time Term from Stabilization");
720
721 timdat.flmpl = inp.GetValue("Lumped Mass Fraction on Left-hand-side");
722 timdat.flmpr = inp.GetValue("Lumped Mass Fraction on Right-hand-side");
723
724 timdat.iCFLworst = 0;
725 if((string)inp.GetValue("Dump CFL") == "True")
726 timdat.iCFLworst = 1;
727
728 intdat.intg[0][0]=inp.GetValue("Quadrature Rule on Interior");
729 intdat.intg[0][1]=inp.GetValue("Quadrature Rule on Boundary");
730 genpar.ibksiz = inp.GetValue("Number of Elements Per Block");
731
732 ((string)inp.GetValue("Turn Off Source Terms for Scalars")
733 == "True") ? sclrs.nosource = 1 : sclrs.nosource = 0;
734 sclrs.tdecay=inp.GetValue("Decay Multiplier for Scalars");
735
736 // TURBULENCE MODELING PARAMETER
737 int tpturb = turbvari.iles-turbvari.irans;
738 int ifrule;
739 if( tpturb != 0 ){
740
741
742 turbvari.nohomog =inp.GetValue("Number of Homogenous Directions");
743
744 if((string)inp.GetValue("Turbulence Wall Model Type") == "Slip Velocity") turbvar.itwmod = 1;
745 else if((string)inp.GetValue("Turbulence Wall Model Type") == "Effective Viscosity") turbvar.itwmod = 2;
746 else turbvar.itwmod = 0;
747 if (turbvari.irans < 0) turbvar.itwmod = turbvar.itwmod*(-1);
748 ifrule = inp.GetValue("Velocity Averaging Steps");
749 turbvar.wtavei =(ifrule >0)? 1.0/ifrule : -1.0/ifrule;
750
751 if(turbvari.iles == 1){
752
753 if((string)inp.GetValue("Dynamic Model Type") == "Bardina") turbvari.iles += 10;
754 else if((string)inp.GetValue("Dynamic Model Type") == "Projection") turbvari.iles += 20;
755
756 ifrule = inp.GetValue("Filter Integration Rule");
757 turbvari.iles += ifrule-1;
758 ifrule = inp.GetValue("Dynamic Model Averaging Steps");
759 turbvar.dtavei = (ifrule >0)? 1.0/ifrule : -1.0/ifrule;
760 turbvar.fwr1 = inp.GetValue("Filter Width Ratio");
761 turbvar.flump = inp.GetValue("Lumping Factor for Filter");
762
763
764 if ((string)inp.GetValue("Model Statistics") == "True" ) {
765 turbvari.modlstats = 1; }
766 else {
767 turbvari.modlstats = 0; }
768
769 if ((string)inp.GetValue("Double Filter") == "True" ) {
770 turbvari.i2filt = 1; }
771 else {
772 turbvari.i2filt = 0; }
773
774 if ((string)inp.GetValue("Model/SUPG Dissipation") == "True" ) {
775 turbvari.idis = 1; }
776 else {
777 turbvari.idis = 0; }
778
779
780 if((string)inp.GetValue("Dynamic Model Type") == "Standard") {
781
782 if((string)inp.GetValue("Dynamic Sub-Model Type") == "None")
783 turbvari.isubmod = 0;
784 else if((string)inp.GetValue("Dynamic Sub-Model Type") =="DFWR")
785 turbvari.isubmod = 1;
786 else if((string)inp.GetValue("Dynamic Sub-Model Type") =="SUPG")
787 turbvari.isubmod = 2;
788 }
789 else if((string)inp.GetValue("Dynamic Model Type") == "Projection") {
790
791 if((string)inp.GetValue("Projection Filter Type") == "Linear")
792 turbvari.ifproj = 0;
793 else if((string)inp.GetValue("Projection Filter Type") =="Quadratic")
794 turbvari.ifproj = 1;
795
796 if((string)inp.GetValue("Dynamic Sub-Model Type") == "None")
797 turbvari.isubmod = 0;
798 else if((string)inp.GetValue("Dynamic Sub-Model Type") =="ConsistentProj")
799 turbvari.isubmod = 1;
800 }
801
802 }
803 }
804
805 // SPEBC MODELING PARAMETERS
806
807 if ( (spebcvr.irscale = inp.GetValue("SPEBC Model Active")) >= 0 ){
808
809 ifrule = inp.GetValue("Velocity Averaging Steps");
810 turbvar.wtavei =(ifrule >0)? 1.0/ifrule : 1.0/inpdat.nstep[0];
811 spebcvr.intpres = inp.GetValue("Interpolate Pressure");
812 spebcvr.plandist = inp.GetValue("Distance between Planes");
813 spebcvr.thetag = inp.GetValue("Theta Angle of Arc");
814 spebcvr.ds = inp.GetValue("Distance for Velocity Averaging");
815 spebcvr.tolerence = inp.GetValue("SPEBC Cylindrical Tolerance");
816 spebcvr.radcyl = inp.GetValue("Radius of recycle plane");
817 spebcvr.rbltin = inp.GetValue("Inlet Boundary Layer Thickness");
818 spebcvr.rvscal = inp.GetValue("Vertical Velocity Scale Factor");
819 }
820
821 // CARDIOVASCULAR MODELING PARAMETERS
822 if ( (string)inp.GetValue("Time Varying Boundary Conditions From File") == "True")
823 nomodule.itvn = 1;
824 else
825 nomodule.itvn = 0;
826
827 if ( nomodule.itvn ==1)
828 nomodule.bcttimescale = inp.GetValue("BCT Time Scale Factor");
829
830 nomodule.ipvsq=0;
831 if( (nomodule.icardio = inp.GetValue("Number of Coupled Surfaces")) ){
832 if ( nomodule.icardio > MAXSURF ) {
833 cout << "Number of Coupled Surfaces > MAXSURF \n";
834 exit(1);
835 }
836 if ( (string)inp.GetValue("Pressure Coupling") == "None")
837 nomodule.ipvsq=0;
838 if ( (string)inp.GetValue("Pressure Coupling") == "Explicit")
839 nomodule.ipvsq=1;
840 if ( (string)inp.GetValue("Pressure Coupling") == "Implicit")
841 nomodule.ipvsq=2;
842 if ( (string)inp.GetValue("Pressure Coupling") == "P-Implicit")
843 nomodule.ipvsq=3;
844
845 if( (nomodule.numResistSrfs=inp.GetValue("Number of Resistance Surfaces")) ){
846 ivec = inp.GetValue("List of Resistance Surfaces");
847 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistResist[i] = 0;
848 for(i=0; i< nomodule.numResistSrfs; i++){
849 nomodule.nsrflistResist[i+1] = ivec[i];
850 }
851 vec = inp.GetValue("Resistance Values");
852 for(i =0; i< MAXSURF+1 ; i++) nomodule.ValueListResist[i] = 0;
853 for(i =0; i< nomodule.numResistSrfs ; i++) nomodule.ValueListResist[i+1] = vec[i];
854 vec.erase(vec.begin(),vec.end());
855 }
856 if( (nomodule.numImpSrfs=inp.GetValue("Number of Impedance Surfaces")) ){
857 ivec = inp.GetValue("List of Impedance Surfaces");
858 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistImp[i] = 0;
859 for(i=0; i< nomodule.numImpSrfs; i++){
860 nomodule.nsrflistImp[i+1] = ivec[i];
861 }
862 if ( (string)inp.GetValue("Impedance From File") == "True")
863 nomodule.impfile = 1; else nomodule.impfile = 0;
864 }
865 if( (nomodule.numRCRSrfs=inp.GetValue("Number of RCR Surfaces")) ){
866 ivec = inp.GetValue("List of RCR Surfaces");
867 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistRCR[i] = 0;
868 for(i=0; i< nomodule.numRCRSrfs; i++){
869 nomodule.nsrflistRCR[i+1] = ivec[i];
870 }
871 if ( (string)inp.GetValue("RCR Values From File") == "True")
872 nomodule.ircrfile = 1; else nomodule.ircrfile = 0;
873 }
874 }
875 nomodule.ideformwall = 0;
876 if((string)inp.GetValue("Deformable Wall")=="True"){
877 nomodule.ideformwall = 1;
878 nomodule.rhovw = inp.GetValue("Density of Vessel Wall");
879 nomodule.thicknessvw = inp.GetValue("Thickness of Vessel Wall");
880 nomodule.evw = inp.GetValue("Young Mod of Vessel Wall");
881 nomodule.rnuvw = inp.GetValue("Poisson Ratio of Vessel Wall");
882 nomodule.rshearconstantvw = inp.GetValue("Shear Constant of Vessel Wall");
883 if((string)inp.GetValue("Wall Mass Matrix for LHS") == "True") nomodule.iwallmassfactor = 1;
884 else nomodule.iwallmassfactor = 0;
885 if((string)inp.GetValue("Wall Stiffness Matrix for LHS") == "True") nomodule.iwallstiffactor = 1;
886 else nomodule.iwallstiffactor = 0;
887 }
888 nomodule.iviscflux = 1;
889 if((string)inp.GetValue("Viscous Flux Flag") == "True") nomodule.iviscflux = 1;
890 if((string)inp.GetValue("Viscous Flux Flag") == "False") nomodule.iviscflux = 0;
891
892
893 // Scaling Parameters Keywords
894
895 outpar.ro = inp.GetValue("Density");
896 outpar.vel = inp.GetValue("Velocity");
897 outpar.press = inp.GetValue("Pressure");
898 outpar.temper = inp.GetValue("Temperature");
899 outpar.entrop = inp.GetValue("Entropy");
900
901 // Step Sequencing
902
903
904 ivec = inp.GetValue("Step Construction");
905 sequence.seqsize = ivec.size();
906 if( sequence.seqsize > 100 || sequence.seqsize < 2 )
907 cerr<<"Sequence size must be between 2 and 100 "<<endl;
908
909 for(i=0; i< sequence.seqsize; i++)
910 sequence.stepseq[i] = ivec[i];
911
912 }
913 catch ( exception &e ) {
914 cout << endl << "Input exception: " << e.what() << endl << endl;
915 ierr = 001;
916 print_error_code(ierr);
917 return ierr;
918 }
919
920 return ierr;
921
922 }
923
print_error_code(int ierr)924 void print_error_code(int ierr) {
925 /*
926 Return Error codes:
927 0xx Input error
928 1xx Solution Control
929 105 Turbulence Model not supported
930
931 2xx Material Properties
932
933 3xx Output Control
934
935 4xx Discretization Control
936
937 5xx Scaling Parameters
938
939 6xx Linear Solver Control
940 601 linear solver type not supported
941 */
942 cout << endl << endl << "Input error detected: " << endl << endl;
943 if ( ierr == 001 ) {
944 cout << endl << "Input Directive not understood" << endl << endl;
945 }
946 if ( ierr == 105 ) {
947 cout << endl << "Turbulence Model Not Supported" << endl << endl;
948 }
949 if ( ierr == 601 ) {
950 cout << endl << "Linear Solver Type Not Supported" << endl << endl;
951 }
952
953 }
954