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