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 if((string)inp.GetValue("Print Statistics") == "True") outpar.ioform = 2; 217 else outpar.ioform = 1; 218 219 if((string)inp.GetValue("Print Wall Fluxes") == "True") outpar.iowflux = 1; 220 else outpar.iowflux = 0; 221 222 if((string)inp.GetValue("Print FieldView") == "True") outpar.iofieldv = 1; 223 else outpar.iofieldv = 0; 224 225 if((string)inp.GetValue("Print ybar") == "True") outpar.ioybar = 1; 226 else outpar.ioybar = 0; 227 228 if((string)inp.GetValue("Print vorticity") == "True") outpar.ivort = 1; 229 else outpar.ivort = 0; 230 231 outpar.nstepsincycle = inp.GetValue("Number of Steps in a Cycle"); 232 outpar.nphasesincycle = inp.GetValue("Number of Phases in a Cycle"); 233 outpar.ncycles_startphaseavg = inp.GetValue("Number of Initial Cycles to Skip in Phase Average"); 234 235 strcpy( outpar.iotype , ((string)inp.GetValue("Data Block Format")).c_str()); 236 strcpy( phasta_iotype , ((string)inp.GetValue("Data Block Format")).c_str()); 237 SONFATH = inp.GetValue("Number of Father Nodes"); 238 239 if((string)inp.GetValue("Print Residual at End of Step") == "True") genpar.lstres = 1; 240 else genpar.lstres = 0; 241 242 if((string)inp.GetValue("Print Error Indicators") == "True") turbvar.ierrcalc = 1; 243 else turbvar.ierrcalc = 0; 244 245 if((string)inp.GetValue("Print Velocity Hessian") == "True") turbvar.ihessian = 1; 246 else turbvar.ihessian = 0; 247 248 if ( turbvar.ierrcalc == 1 ) 249 turbvari.ierrsmooth = inp.GetValue("Number of Error Smoothing Iterations"); 250 251 for(i=0;i<MAXSURF+1; i++) aerfrc.nsrflist[i] = 0; 252 int nsrfCM = inp.GetValue("Number of Force Surfaces"); 253 if (nsrfCM > 0) { 254 vector<int> ivec = inp.GetValue("Surface ID's for Force Calculation"); 255 for(i=0; i< nsrfCM; i++){ 256 aerfrc.nsrflist[ivec[i]] = 1; 257 // cout <<"surface in force list "<< ivec[i] << endl; 258 } 259 ivec.erase(ivec.begin(),ivec.end()); 260 } 261 262 aerfrc.isrfIM = inp.GetValue("Surface ID for Integrated Mass"); 263 264 outpar.iv_rankpercore = inp.GetValue("Ranks per core"); 265 outpar.iv_corepernode = inp.GetValue("Cores per node"); 266 267 turbvari.iramp=0; 268 if((string)inp.GetValue("Ramp Inflow") == "True") turbvari.iramp=1; 269 if(turbvari.iramp == 1) { 270 vec = inp.GetValue("Mdot Ramp Inflow Start and Stop"); 271 for(i=0; i<2 ; i++){ 272 turbvar.rampmdot[0][i]=vec[i]; 273 } 274 vec = inp.GetValue("Mdot Ramp Lower FC Start and Stop"); 275 for(i=0; i<2 ; i++){ 276 turbvar.rampmdot[1][i]=vec[i]; 277 } 278 vec = inp.GetValue("Mdot Ramp Upper FC Start and Stop"); 279 for(i=0; i<2 ; i++){ 280 turbvar.rampmdot[2][i]=vec[i]; 281 } 282 } 283 284 //Limiting 285 vec = inp.GetValue("Limit u1"); 286 for(i=0; i<3 ; i++){ 287 turbvar.ylimit[0][i] = vec[i]; 288 } 289 vec.erase(vec.begin(),vec.end()); 290 291 vec = inp.GetValue("Limit u2"); 292 for(i=0; i<3 ; i++){ 293 turbvar.ylimit[1][i] = vec[i]; 294 } 295 vec.erase(vec.begin(),vec.end()); 296 297 vec = inp.GetValue("Limit u3"); 298 for(i=0; i<3 ; i++){ 299 turbvar.ylimit[2][i] = vec[i]; 300 } 301 vec.erase(vec.begin(),vec.end()); 302 303 vec = inp.GetValue("Limit Pressure"); 304 for(i=0; i<3 ; i++){ 305 turbvar.ylimit[3][i] = vec[i]; 306 } 307 vec.erase(vec.begin(),vec.end()); 308 309 vec = inp.GetValue("Limit Temperature"); 310 for(i=0; i<3 ; i++){ 311 turbvar.ylimit[4][i] = vec[i]; 312 } 313 vec.erase(vec.begin(),vec.end()); 314 315 //Material Properties Keywords 316 matdat.nummat = levlset.iLSet+1; 317 if((string)inp.GetValue("Shear Law") == "Constant Viscosity") 318 for(i=0; i < levlset.iLSet+1; i++) matdat.matflg[i][1] = 0; 319 320 if((string)inp.GetValue("Bulk Viscosity Law") == "Constant Bulk Viscosity") 321 for(i=0; i < levlset.iLSet+1; i++) matdat.matflg[i][2] = 0; 322 323 mmatpar.pr = inp.GetValue("Prandtl Number"); 324 325 if((string)inp.GetValue("Conductivity Law") == "Constant Conductivity") 326 for(i=0; i < levlset.iLSet+1; i++) matdat.matflg[i][3] = 0; 327 328 vec = inp.GetValue("Density"); 329 for(i=0; i< levlset.iLSet +1 ; i++){ 330 matdat.datmat[i][0][0] = vec[i]; 331 } 332 vec.erase(vec.begin(),vec.end()); 333 334 vec = inp.GetValue("Viscosity"); 335 for(i=0; i< levlset.iLSet +1 ; i++){ 336 matdat.datmat[i][1][0] = vec[i]; 337 } 338 vec.erase(vec.begin(),vec.end()); 339 340 // vec = inp.GetValue("Specific Heat"); 341 for(i=0; i< levlset.iLSet +1 ; i++){ 342 matdat.datmat[i][2][0] = 0; 343 } 344 // vec.erase(vec.begin(),vec.end()); 345 346 vec = inp.GetValue("Thermal Conductivity"); 347 for(i=0; i< levlset.iLSet +1 ; i++){ 348 matdat.datmat[i][3][0] = vec[i]; 349 } 350 vec.erase(vec.begin(),vec.end()); 351 352 vec = inp.GetValue("Scalar Diffusivity"); 353 for(i=0; i< solscalr ; i++){ 354 sclrs.scdiff[i] = vec[i]; 355 } 356 vec.erase(vec.begin(),vec.end()); 357 358 if((string)inp.GetValue("Zero Mean Pressure") == "True") 359 turbvar.pzero=1; 360 361 turbvar.rmutarget = inp.GetValue("Target Viscosity For Step NSTEP"); 362 363 if ( (string)inp.GetValue("Body Force Option") == "None" ) { 364 for( i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 0; 365 } 366 else if ( (string)inp.GetValue("Body Force Option") == "Vector" ) { 367 for( i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 1; 368 } 369 else if ( (string)inp.GetValue("Body Force Option") == "User e3source.f" ) { 370 for( i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 3; 371 } 372 else if ( (string)inp.GetValue("Body Force Option") == "Boussinesq" ) { 373 for(i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 2; 374 } 375 else if ( (string)inp.GetValue("Body Force Option") == "Cooling Analytic" ) { 376 for(i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 4; 377 } 378 else if ( (string)inp.GetValue("Body Force Option") == "Cooling Initial Condition" ) { 379 for(i=0; i< levlset.iLSet +1 ; i++) matdat.matflg[i][4] = 5; 380 } 381 382 // the following block of stuff is common to all cooling type sponges. 383 // Specific stuff belongs in the conditionals above 384 385 if(matdat.matflg[0][4] >=4) { 386 spongevar.betamax = inp.GetValue("Maximum Value of Sponge Parameter"); 387 spongevar.zinsponge = inp.GetValue("Inflow Cooling Sponge Ends at z"); 388 spongevar.zoutsponge= inp.GetValue("Outflow Cooling Sponge Begins at z"); 389 spongevar.radsponge = inp.GetValue("Radial Cooling Sponge Begins at r"); 390 spongevar.grthosponge = inp.GetValue("Sponge Growth Coefficient Outflow"); 391 spongevar.grthisponge = inp.GetValue("Sponge Growth Coefficient Inflow"); 392 393 394 spongevar.spongecontinuity = 0; 395 spongevar.spongemomentum1 = 0; 396 spongevar.spongemomentum2 = 0; 397 spongevar.spongemomentum3 = 0; 398 spongevar.spongeenergy = 0; 399 400 if((string)inp.GetValue("Sponge for Continuity Equation") == "True") 401 spongevar.spongecontinuity = 1; 402 if((string)inp.GetValue("Sponge for x Momentum Equation") == "True") 403 spongevar.spongemomentum1 = 1; 404 if((string)inp.GetValue("Sponge for y Momentum Equation") == "True") 405 spongevar.spongemomentum2 = 1; 406 if((string)inp.GetValue("Sponge for z Momentum Equation") == "True") 407 spongevar.spongemomentum3 = 1; 408 if((string)inp.GetValue("Sponge for Energy Equation") == "True") 409 spongevar.spongeenergy = 1; 410 411 } 412 413 vec = inp.GetValue("Body Force"); 414 for(i=0; i< levlset.iLSet +1 ; i++){ 415 matdat.datmat[i][4][0] = vec[0+i*3]; 416 matdat.datmat[i][4][1] = vec[1+i*3]; 417 matdat.datmat[i][4][2] = vec[2+i*3]; 418 } 419 vec.erase(vec.begin(),vec.end()); 420 421 vec = inp.GetValue("Body Force Pressure Gradient"); 422 for(i=0; i< levlset.iLSet +1 ; i++){ 423 matdat.datmat[i][6][0] = vec[0+i*3]; 424 matdat.datmat[i][6][1] = vec[1+i*3]; 425 matdat.datmat[i][6][2] = vec[2+i*3]; 426 } 427 vec.erase(vec.begin(),vec.end()); 428 429 if ( (string)inp.GetValue("Surface Tension Option") == "No" ){ 430 genpar.isurf = 0; 431 } 432 else if ((string)inp.GetValue("Surface Tension Option") == "Yes" ){ 433 genpar.isurf = 1; 434 } 435 else { 436 cout << " Surface Tension: Only Legal Values (Yes, No) "; 437 cout << endl; 438 exit(1); 439 } 440 if( genpar.isurf > 0) { 441 genpar.Bo = inp.GetValue("Bond Number"); 442 } 443 444 genpar.EntropyPressure = inp.GetValue("Entropy Form of Pressure Constraint on Weight Space"); 445 446 447 if ( (string)inp.GetValue("Rotating Frame of Reference") == "True" ) { 448 matdat.matflg[0][5] = 1; 449 vec = inp.GetValue("Rotating Frame of Reference Rotation Rate"); 450 matdat.datmat[0][5][0] = vec[0]; 451 matdat.datmat[0][5][1] = vec[1]; 452 matdat.datmat[0][5][2] = vec[2]; 453 vec.erase(vec.begin(),vec.end()); 454 } 455 else { 456 matdat.matflg[0][5] = 0; 457 matdat.datmat[0][5][0] = 0.; 458 matdat.datmat[0][5][1] = 0.; 459 matdat.datmat[0][5][2] = 0.; 460 } 461 462 463 //Linear Solver parameters 464 if( (string)inp.GetValue("Solver Type") =="ACUSIM with P Projection" ){ 465 incomp.iprjFlag = 0; incomp.ipresPrjFlag=1;} 466 else if ( (string)inp.GetValue("Solver Type") =="ACUSIM" ){ 467 incomp.iprjFlag = 0; incomp.ipresPrjFlag=0;} 468 else if( (string)inp.GetValue("Solver Type") =="ACUSIM with Velocity Projection" ){ 469 incomp.iprjFlag = 1; incomp.ipresPrjFlag=0;} 470 else if( (string)inp.GetValue("Solver Type") =="ACUSIM with Full Projection" ){ 471 incomp.iprjFlag = 1; incomp.ipresPrjFlag=1;} 472 else if( (string)inp.GetValue("Solver Type") =="GMRES Matrix Free"){ 473 inpdat.impl[0] += 10*solflow;} 474 else if( (string)inp.GetValue("Solver Type") =="GMRES EBE"){ 475 inpdat.impl[0] += 20*solflow;} 476 //GMRES sparse is assumed default and has the value of 10, MFG 20, 477 // EBE 30 478 479 480 // inpdat.niter[0] = inp.GetValue("Number of Solves per Time Step"); 481 solpar.nGMRES = inp.GetValue("Number of GMRES Sweeps per Solve"); 482 solpar.Kspace = inp.GetValue("Number of Krylov Vectors per GMRES Sweep"); 483 inpdat.LHSupd[0] = inp.GetValue("Number of Solves per Left-hand-side Formation"); 484 inpdat.epstol[0] = inp.GetValue("Tolerance on Momentum Equations"); 485 incomp.prestol = inp.GetValue("Tolerance on ACUSIM Pressure Projection"); 486 incomp.minIters = inp.GetValue("Minimum Number of Iterations per Nonlinear Iteration"); 487 incomp.maxIters = inp.GetValue("Maximum Number of Iterations per Nonlinear Iteration"); 488 inpdat.deltol[0][0]=inp.GetValue("Velocity Delta Ratio"); 489 inpdat.deltol[1][0]=inp.GetValue("Pressure Delta Ratio"); 490 incomp.nPrjs = inp.GetValue("Number of Velocity Projection Vectors"); 491 incomp.nPresPrjs = inp.GetValue("Number of Pressure Projection Vectors"); 492 incomp.iverbose = inp.GetValue("ACUSIM Verbosity Level"); 493 494 if(solheat==1){ 495 inpdat.epstol[1]=inp.GetValue("Temperature Solver Tolerance"); 496 inpdat.LHSupd[1]=inp.GetValue("Number of Solves of Temperature per Left-hand-side Formation"); 497 } 498 499 // The following is where you should put any inputs that are able to 500 // input differently for each scalar. It is a little tedious in the code 501 // but it should make the solver.inp easier to understand. Note this will 502 // require some care with regression tests. 503 504 505 if(solscalr>0){ 506 inpdat.epstol[2]=inp.GetValue("Scalar 1 Solver Tolerance"); 507 inpdat.LHSupd[2]=inp.GetValue("Number of Solves of Scalar 1 per Left-hand-side Formation"); 508 509 vec = inp.GetValue("Limit Scalar 1"); 510 for(i=0; i<3 ; i++){ 511 turbvar.ylimit[5][i] = vec[i]; 512 } 513 vec.erase(vec.begin(),vec.end()); 514 } 515 516 if(solscalr>1){ 517 inpdat.epstol[3]=inp.GetValue("Scalar 2 Solver Tolerance"); 518 inpdat.LHSupd[3]=inp.GetValue("Number of Solves of Scalar 2 per Left-hand-side Formation"); 519 520 vec = inp.GetValue("Limit Scalar 2"); 521 for(i=0; i<3 ; i++){ 522 turbvar.ylimit[6][i] = vec[i]; 523 } 524 vec.erase(vec.begin(),vec.end()); 525 } 526 527 if(solscalr>2){ 528 inpdat.epstol[4]=inp.GetValue("Scalar 3 Solver Tolerance"); 529 inpdat.LHSupd[4]=inp.GetValue("Number of Solves of Scalar 3 per Left-hand-side Formation"); 530 531 vec = inp.GetValue("Limit Scalar 3"); 532 for(i=0; i<3 ; i++){ 533 turbvar.ylimit[7][i] = vec[i]; 534 } 535 vec.erase(vec.begin(),vec.end()); 536 } 537 538 if(solscalr>3){ 539 inpdat.epstol[5]=inp.GetValue("Scalar 4 Solver Tolerance"); 540 inpdat.LHSupd[5]=inp.GetValue("Number of Solves of Scalar 4 per Left-hand-side Formation"); 541 542 vec = inp.GetValue("Limit Scalar 4"); 543 for(i=0; i<3 ; i++){ 544 turbvar.ylimit[8][i] = vec[i]; 545 } 546 vec.erase(vec.begin(),vec.end()); 547 } 548 549 // DISCRETIZATION CONTROL 550 551 genpar.ipord = inp.GetValue("Basis Function Order"); 552 if((string)inp.GetValue("Time Integration Rule") == "First Order") 553 inpdat.rhoinf[0] = -1 ; 554 else inpdat.rhoinf[0] = (double)inp.GetValue("Time Integration Rho Infinity"); 555 if((string)inp.GetValue("Predictor at Start of Step")=="Same Velocity") 556 genpar.ipred = 1; 557 if((string)inp.GetValue("Predictor at Start of Step")=="Zero Acceleration") 558 genpar.ipred = 2; 559 if((string)inp.GetValue("Predictor at Start of Step")=="Same Acceleration") 560 genpar.ipred = 3; 561 if((string)inp.GetValue("Predictor at Start of Step")=="Same Delta") 562 genpar.ipred = 4; 563 564 if((string)inp.GetValue("Weak Form") == "Galerkin") 565 solpar.ivart = 1; 566 if((string)inp.GetValue("Weak Form") == "SUPG") 567 solpar.ivart = 2; 568 569 if((string)inp.GetValue("Flow Advection Form") == "Convective") 570 solpar.iconvflow = 2; 571 else if((string)inp.GetValue("Flow Advection Form") == "Conservative") 572 solpar.iconvflow = 1; 573 if((string)inp.GetValue("Scalar Advection Form") == "Convective") 574 solpar.iconvsclr = 2; 575 else if((string)inp.GetValue("Scalar Advection Form") == "Conservative") 576 solpar.iconvsclr = 1; 577 if((string)inp.GetValue("Use Conservative Scalar Convection Velocity") == "True") 578 sclrs.consrv_sclr_conv_vel = 1; 579 else if((string)inp.GetValue("Use Conservative Scalar Convection Velocity") == "False") 580 sclrs.consrv_sclr_conv_vel = 0; 581 // TAU INPUT 582 if((string)inp.GetValue("Tau Matrix") == "Diagonal-Shakib") 583 genpar.itau = 0; 584 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Franca") 585 genpar.itau =1; 586 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Jansen(dev)") 587 genpar.itau = 2; 588 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Compressible") 589 genpar.itau = 3; 590 else if((string)inp.GetValue("Tau Matrix") == "Matrix-Mallet") 591 genpar.itau = 10; 592 else if((string)inp.GetValue("Tau Matrix") == "Matrix-Modal") 593 genpar.itau = 11; 594 595 genpar.dtsfct = inp.GetValue("Tau Time Constant"); 596 genpar.taucfct = inp.GetValue("Tau C Scale Factor"); 597 598 // FLOW DISCONTINUITY CAPTURING 599 600 if((string)inp.GetValue("Discontinuity Capturing") == "Off") solpar.iDC = 0; 601 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-mallet") solpar.iDC = 1; 602 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-quadratic") solpar.iDC = 2; 603 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-minimum") solpar.iDC = 3; 604 else { 605 cout<< "Condition not defined for Discontinuity Capturing \n "; 606 exit(1); 607 } 608 609 // SCALAR DISCONTINUITY CAPTURING 610 611 vector<int> ivec = inp.GetValue("Scalar Discontinuity Capturing"); 612 for(i=0; i< 2; i++) solpar.idcsclr[i] = ivec[i]; 613 ivec.erase(ivec.begin(),ivec.end()); 614 615 616 // if((string)inp.GetValue("Scalar Discontinuity Capturing") == "No") solpar.idcsclr = 0; 617 // else if((string)inp.GetValue("Scalar Discontinuity Capturing") == "1") solpar.idcsclr = 1; 618 // else if((string)inp.GetValue("Scalar Discontinuity Capturing") == "2") solpar.idcsclr = 2; 619 // else { 620 // cout<< "Condition not defined for Scalar Discontinuity Capturing \n "; 621 // exit(1); 622 // } 623 if((string)inp.GetValue("Include Viscous Correction in Stabilization") == "True") 624 { 625 if(genpar.ipord == 1 ) genpar.idiff = 1; 626 else genpar.idiff = 2; 627 } 628 else { genpar.idiff = 0;} 629 630 genpar.irampViscOutlet = (int)inp.GetValue("Ramp Up Viscosity Near Outlet"); 631 632 genpar.istretchOutlet = (int)inp.GetValue("Stretch X Coordinate Near Outlet"); 633 634 genpar.iremoveStabTimeTerm = (int)inp.GetValue("Remove Time Term from Stabilization"); 635 636 timdat.flmpl = inp.GetValue("Lumped Mass Fraction on Left-hand-side"); 637 timdat.flmpr = inp.GetValue("Lumped Mass Fraction on Right-hand-side"); 638 639 timdat.iCFLworst = 0; 640 if((string)inp.GetValue("Dump CFL") == "True") 641 timdat.iCFLworst = 1; 642 643 intdat.intg[0][0]=inp.GetValue("Quadrature Rule on Interior"); 644 intdat.intg[0][1]=inp.GetValue("Quadrature Rule on Boundary"); 645 genpar.ibksiz = inp.GetValue("Number of Elements Per Block"); 646 647 ((string)inp.GetValue("Turn Off Source Terms for Scalars") 648 == "True")? sclrs.nosource=1:sclrs.nosource=0; 649 sclrs.tdecay=inp.GetValue("Decay Multiplier for Scalars"); 650 651 // TURBULENCE MODELING PARAMETER 652 int tpturb = turbvari.iles-turbvari.irans; 653 int ifrule; 654 if( tpturb != 0 ){ 655 656 657 turbvari.nohomog =inp.GetValue("Number of Homogenous Directions"); 658 659 if((string)inp.GetValue("Turbulence Wall Model Type") == "Slip Velocity") turbvar.itwmod = 1; 660 else if((string)inp.GetValue("Turbulence Wall Model Type") == "Effective Viscosity") turbvar.itwmod = 2; 661 else turbvar.itwmod = 0; 662 if (turbvari.irans < 0) turbvar.itwmod = turbvar.itwmod*(-1); 663 ifrule = inp.GetValue("Velocity Averaging Steps"); 664 turbvar.wtavei =(ifrule >0)? 1.0/ifrule : -1.0/ifrule; 665 666 if(turbvari.iles == 1){ 667 668 if((string)inp.GetValue("Dynamic Model Type") == "Bardina") turbvari.iles += 10; 669 else if((string)inp.GetValue("Dynamic Model Type") == "Projection") turbvari.iles += 20; 670 671 ifrule = inp.GetValue("Filter Integration Rule"); 672 turbvari.iles += ifrule-1; 673 ifrule = inp.GetValue("Dynamic Model Averaging Steps"); 674 turbvar.dtavei = (ifrule >0)? 1.0/ifrule : -1.0/ifrule; 675 turbvar.fwr1 = inp.GetValue("Filter Width Ratio"); 676 turbvar.flump = inp.GetValue("Lumping Factor for Filter"); 677 678 679 if ((string)inp.GetValue("Model Statistics") == "True" ) { 680 turbvari.modlstats = 1; } 681 else { 682 turbvari.modlstats = 0; } 683 684 if ((string)inp.GetValue("Double Filter") == "True" ) { 685 turbvari.i2filt = 1; } 686 else { 687 turbvari.i2filt = 0; } 688 689 if ((string)inp.GetValue("Model/SUPG Dissipation") == "True" ) { 690 turbvari.idis = 1; } 691 else { 692 turbvari.idis = 0; } 693 694 695 if((string)inp.GetValue("Dynamic Model Type") == "Standard") { 696 697 if((string)inp.GetValue("Dynamic Sub-Model Type") == "None") 698 turbvari.isubmod = 0; 699 else if((string)inp.GetValue("Dynamic Sub-Model Type") =="DFWR") 700 turbvari.isubmod = 1; 701 else if((string)inp.GetValue("Dynamic Sub-Model Type") =="SUPG") 702 turbvari.isubmod = 2; 703 } 704 else if((string)inp.GetValue("Dynamic Model Type") == "Projection") { 705 706 if((string)inp.GetValue("Projection Filter Type") == "Linear") 707 turbvari.ifproj = 0; 708 else if((string)inp.GetValue("Projection Filter Type") =="Quadratic") 709 turbvari.ifproj = 1; 710 711 if((string)inp.GetValue("Dynamic Sub-Model Type") == "None") 712 turbvari.isubmod = 0; 713 else if((string)inp.GetValue("Dynamic Sub-Model Type") =="ConsistentProj") 714 turbvari.isubmod = 1; 715 } 716 717 } 718 } 719 720 // SPEBC MODELING PARAMETERS 721 722 if ( (spebcvr.irscale = inp.GetValue("SPEBC Model Active")) >= 0 ){ 723 724 ifrule = inp.GetValue("Velocity Averaging Steps"); 725 turbvar.wtavei =(ifrule >0)? 1.0/ifrule : 1.0/inpdat.nstep[0]; 726 spebcvr.intpres = inp.GetValue("Interpolate Pressure"); 727 spebcvr.plandist = inp.GetValue("Distance between Planes"); 728 spebcvr.thetag = inp.GetValue("Theta Angle of Arc"); 729 spebcvr.ds = inp.GetValue("Distance for Velocity Averaging"); 730 spebcvr.tolerence = inp.GetValue("SPEBC Cylindrical Tolerance"); 731 spebcvr.radcyl = inp.GetValue("Radius of recycle plane"); 732 spebcvr.rbltin = inp.GetValue("Inlet Boundary Layer Thickness"); 733 spebcvr.rvscal = inp.GetValue("Vertical Velocity Scale Factor"); 734 } 735 736 // CARDIOVASCULAR MODELING PARAMETERS 737 if ( (string)inp.GetValue("Time Varying Boundary Conditions From File") == "True") 738 nomodule.itvn = 1; 739 else 740 nomodule.itvn = 0; 741 742 if ( nomodule.itvn ==1) 743 nomodule.bcttimescale = inp.GetValue("BCT Time Scale Factor"); 744 745 nomodule.ipvsq=0; 746 if(nomodule.icardio = inp.GetValue("Number of Coupled Surfaces")){ 747 if ( nomodule.icardio > MAXSURF ) { 748 cout << "Number of Coupled Surfaces > MAXSURF \n"; 749 exit(1); 750 } 751 if ( (string)inp.GetValue("Pressure Coupling") == "None") 752 nomodule.ipvsq=0; 753 if ( (string)inp.GetValue("Pressure Coupling") == "Explicit") 754 nomodule.ipvsq=1; 755 if ( (string)inp.GetValue("Pressure Coupling") == "Implicit") 756 nomodule.ipvsq=2; 757 if ( (string)inp.GetValue("Pressure Coupling") == "P-Implicit") 758 nomodule.ipvsq=3; 759 760 if(nomodule.numResistSrfs=inp.GetValue("Number of Resistance Surfaces")){ 761 ivec = inp.GetValue("List of Resistance Surfaces"); 762 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistResist[i] = 0; 763 for(i=0; i< nomodule.numResistSrfs; i++){ 764 nomodule.nsrflistResist[i+1] = ivec[i]; 765 } 766 vec = inp.GetValue("Resistance Values"); 767 for(i =0; i< MAXSURF+1 ; i++) nomodule.ValueListResist[i] = 0; 768 for(i =0; i< nomodule.numResistSrfs ; i++) nomodule.ValueListResist[i+1] = vec[i]; 769 vec.erase(vec.begin(),vec.end()); 770 } 771 if(nomodule.numImpSrfs=inp.GetValue("Number of Impedance Surfaces")){ 772 ivec = inp.GetValue("List of Impedance Surfaces"); 773 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistImp[i] = 0; 774 for(i=0; i< nomodule.numImpSrfs; i++){ 775 nomodule.nsrflistImp[i+1] = ivec[i]; 776 } 777 if ( (string)inp.GetValue("Impedance From File") == "True") 778 nomodule.impfile = 1; else nomodule.impfile = 0; 779 } 780 if(nomodule.numRCRSrfs=inp.GetValue("Number of RCR Surfaces")){ 781 ivec = inp.GetValue("List of RCR Surfaces"); 782 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistRCR[i] = 0; 783 for(i=0; i< nomodule.numRCRSrfs; i++){ 784 nomodule.nsrflistRCR[i+1] = ivec[i]; 785 } 786 if ( (string)inp.GetValue("RCR Values From File") == "True") 787 nomodule.ircrfile = 1; else nomodule.ircrfile = 0; 788 } 789 } 790 nomodule.ideformwall = 0; 791 if((string)inp.GetValue("Deformable Wall")=="True"){ 792 nomodule.ideformwall = 1; 793 nomodule.rhovw = inp.GetValue("Density of Vessel Wall"); 794 nomodule.thicknessvw = inp.GetValue("Thickness of Vessel Wall"); 795 nomodule.evw = inp.GetValue("Young Mod of Vessel Wall"); 796 nomodule.rnuvw = inp.GetValue("Poisson Ratio of Vessel Wall"); 797 nomodule.rshearconstantvw = inp.GetValue("Shear Constant of Vessel Wall"); 798 if((string)inp.GetValue("Wall Mass Matrix for LHS") == "True") nomodule.iwallmassfactor = 1; 799 else nomodule.iwallmassfactor = 0; 800 if((string)inp.GetValue("Wall Stiffness Matrix for LHS") == "True") nomodule.iwallstiffactor = 1; 801 else nomodule.iwallstiffactor = 0; 802 } 803 nomodule.iviscflux = 1; 804 if((string)inp.GetValue("Viscous Flux Flag") == "True") nomodule.iviscflux = 1; 805 if((string)inp.GetValue("Viscous Flux Flag") == "False") nomodule.iviscflux = 0; 806 807 808 // Scaling Parameters Keywords 809 810 outpar.ro = inp.GetValue("Density"); 811 outpar.vel = inp.GetValue("Velocity"); 812 outpar.press = inp.GetValue("Pressure"); 813 outpar.temper = inp.GetValue("Temperature"); 814 outpar.entrop = inp.GetValue("Entropy"); 815 816 // Step Sequencing 817 818 819 ivec = inp.GetValue("Step Construction"); 820 sequence.seqsize = ivec.size(); 821 if( sequence.seqsize > 100 || sequence.seqsize < 2 ) 822 cerr<<"Sequence size must be between 2 and 100 "<<endl; 823 824 for(i=0; i< sequence.seqsize; i++) 825 sequence.stepseq[i] = ivec[i]; 826 827 } 828 catch ( exception &e ) { 829 cout << endl << "Input exception: " << e.what() << endl << endl; 830 ierr = 001; 831 print_error_code(ierr); 832 return ierr; 833 } 834 835 return ierr; 836 837 } 838 839 void print_error_code(int ierr) { 840 /* 841 Return Error codes: 842 0xx Input error 843 1xx Solution Control 844 105 Turbulence Model not supported 845 846 2xx Material Properties 847 848 3xx Output Control 849 850 4xx Discretization Control 851 852 5xx Scaling Parameters 853 854 6xx Linear Solver Control 855 601 linear solver type not supported 856 */ 857 cout << endl << endl << "Input error detected: " << endl << endl; 858 if ( ierr == 001 ) { 859 cout << endl << "Input Directive not understood" << endl << endl; 860 } 861 if ( ierr == 105 ) { 862 cout << endl << "Turbulence Model Not Supported" << endl << endl; 863 } 864 if ( ierr == 601 ) { 865 cout << endl << "Linear Solver Type Not Supported" << endl << endl; 866 } 867 868 } 869