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