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 incomp.iprjFlag = 0; incomp.ipresPrjFlag=0; inpdat.svLSFlag=0; 448 if( (string)inp.GetValue("Solver Type") =="svLS" ){ 449 inpdat.svLSFlag = 1; } 450 if( (string)inp.GetValue("Solver Type") =="ACUSIM with P Projection" ){ 451 incomp.iprjFlag = 0; incomp.ipresPrjFlag=1;} 452 else if ( (string)inp.GetValue("Solver Type") =="ACUSIM" ){ 453 incomp.iprjFlag = 0; incomp.ipresPrjFlag=0;} 454 else if( (string)inp.GetValue("Solver Type") =="ACUSIM with Velocity Projection" ){ 455 incomp.iprjFlag = 1; incomp.ipresPrjFlag=0;} 456 else if( (string)inp.GetValue("Solver Type") =="ACUSIM with Full Projection" ){ 457 incomp.iprjFlag = 1; incomp.ipresPrjFlag=1;} 458 else if( (string)inp.GetValue("Solver Type") =="GMRES Matrix Free"){ 459 inpdat.impl[0] += 10*solflow;} 460 else if( (string)inp.GetValue("Solver Type") =="GMRES EBE"){ 461 inpdat.impl[0] += 20*solflow;} 462 //GMRES sparse is assumed default and has the value of 10, MFG 20, 463 // EBE 30 464 465 466 // inpdat.niter[0] = inp.GetValue("Number of Solves per Time Step"); 467 solpar.nGMRES = inp.GetValue("Number of GMRES Sweeps per Solve"); 468 solpar.Kspace = inp.GetValue("Number of Krylov Vectors per GMRES Sweep"); 469 inpdat.LHSupd[0] = inp.GetValue("Number of Solves per Left-hand-side Formation"); 470 inpdat.epstol[0] = inp.GetValue("Tolerance on Momentum Equations"); 471 incomp.prestol = inp.GetValue("Tolerance on ACUSIM Pressure Projection"); 472 incomp.minIters = inp.GetValue("Minimum Number of Iterations per Nonlinear Iteration"); 473 incomp.maxIters = inp.GetValue("Maximum Number of Iterations per Nonlinear Iteration"); 474 inpdat.deltol[0][0]=inp.GetValue("Velocity Delta Ratio"); 475 inpdat.deltol[1][0]=inp.GetValue("Pressure Delta Ratio"); 476 incomp.nPrjs = inp.GetValue("Number of Velocity Projection Vectors"); 477 incomp.nPresPrjs = inp.GetValue("Number of Pressure Projection Vectors"); 478 incomp.iverbose = inp.GetValue("ACUSIM Verbosity Level"); 479 480 if(solheat==1){ 481 inpdat.epstol[1]=inp.GetValue("Temperature Solver Tolerance"); 482 inpdat.LHSupd[1]=inp.GetValue("Number of Solves of Temperature per Left-hand-side Formation"); 483 } 484 485 // The following is where you should put any inputs that are able to 486 // input differently for each scalar. It is a little tedious in the code 487 // but it should make the solver.inp easier to understand. Note this will 488 // require some care with regression tests. 489 490 491 if(solscalr>0){ 492 inpdat.epstol[2]=inp.GetValue("Scalar 1 Solver Tolerance"); 493 inpdat.LHSupd[2]=inp.GetValue("Number of Solves of Scalar 1 per Left-hand-side Formation"); 494 495 vec = inp.GetValue("Limit Scalar 1"); 496 for(i=0; i<3 ; i++){ 497 turbvar.ylimit[5][i] = vec[i]; 498 } 499 vec.erase(vec.begin(),vec.end()); 500 } 501 502 if(solscalr>1){ 503 inpdat.epstol[3]=inp.GetValue("Scalar 2 Solver Tolerance"); 504 inpdat.LHSupd[3]=inp.GetValue("Number of Solves of Scalar 2 per Left-hand-side Formation"); 505 506 vec = inp.GetValue("Limit Scalar 2"); 507 for(i=0; i<3 ; i++){ 508 turbvar.ylimit[6][i] = vec[i]; 509 } 510 vec.erase(vec.begin(),vec.end()); 511 } 512 513 if(solscalr>2){ 514 inpdat.epstol[4]=inp.GetValue("Scalar 3 Solver Tolerance"); 515 inpdat.LHSupd[4]=inp.GetValue("Number of Solves of Scalar 3 per Left-hand-side Formation"); 516 517 vec = inp.GetValue("Limit Scalar 3"); 518 for(i=0; i<3 ; i++){ 519 turbvar.ylimit[7][i] = vec[i]; 520 } 521 vec.erase(vec.begin(),vec.end()); 522 } 523 524 if(solscalr>3){ 525 inpdat.epstol[5]=inp.GetValue("Scalar 4 Solver Tolerance"); 526 inpdat.LHSupd[5]=inp.GetValue("Number of Solves of Scalar 4 per Left-hand-side Formation"); 527 528 vec = inp.GetValue("Limit Scalar 4"); 529 for(i=0; i<3 ; i++){ 530 turbvar.ylimit[8][i] = vec[i]; 531 } 532 vec.erase(vec.begin(),vec.end()); 533 } 534 535 // DISCRETIZATION CONTROL 536 537 genpar.ipord = inp.GetValue("Basis Function Order"); 538 if((string)inp.GetValue("Time Integration Rule") == "First Order") 539 inpdat.rhoinf[0] = -1 ; 540 else inpdat.rhoinf[0] = (double)inp.GetValue("Time Integration Rho Infinity"); 541 if((string)inp.GetValue("Predictor at Start of Step")=="Same Velocity") 542 genpar.ipred = 1; 543 if((string)inp.GetValue("Predictor at Start of Step")=="Zero Acceleration") 544 genpar.ipred = 2; 545 if((string)inp.GetValue("Predictor at Start of Step")=="Same Acceleration") 546 genpar.ipred = 3; 547 if((string)inp.GetValue("Predictor at Start of Step")=="Same Delta") 548 genpar.ipred = 4; 549 550 if((string)inp.GetValue("Weak Form") == "Galerkin") 551 solpar.ivart = 1; 552 if((string)inp.GetValue("Weak Form") == "SUPG") 553 solpar.ivart = 2; 554 555 if((string)inp.GetValue("Flow Advection Form") == "Convective") 556 solpar.iconvflow = 2; 557 else if((string)inp.GetValue("Flow Advection Form") == "Conservative") 558 solpar.iconvflow = 1; 559 if((string)inp.GetValue("Scalar Advection Form") == "Convective") 560 solpar.iconvsclr = 2; 561 else if((string)inp.GetValue("Scalar Advection Form") == "Conservative") 562 solpar.iconvsclr = 1; 563 if((string)inp.GetValue("Use Conservative Scalar Convection Velocity") == "True") 564 sclrs.consrv_sclr_conv_vel = 1; 565 else if((string)inp.GetValue("Use Conservative Scalar Convection Velocity") == "False") 566 sclrs.consrv_sclr_conv_vel = 0; 567 // TAU INPUT 568 if((string)inp.GetValue("Tau Matrix") == "Diagonal-Shakib") 569 genpar.itau = 0; 570 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Franca") 571 genpar.itau =1; 572 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Jansen(dev)") 573 genpar.itau = 2; 574 else if((string)inp.GetValue("Tau Matrix") == "Diagonal-Compressible") 575 genpar.itau = 3; 576 else if((string)inp.GetValue("Tau Matrix") == "Matrix-Mallet") 577 genpar.itau = 10; 578 else if((string)inp.GetValue("Tau Matrix") == "Matrix-Modal") 579 genpar.itau = 11; 580 581 genpar.dtsfct = inp.GetValue("Tau Time Constant"); 582 genpar.taucfct = inp.GetValue("Tau C Scale Factor"); 583 584 // FLOW DISCONTINUITY CAPTURING 585 586 if((string)inp.GetValue("Discontinuity Capturing") == "Off") solpar.iDC = 0; 587 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-mallet") solpar.iDC = 1; 588 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-quadratic") solpar.iDC = 2; 589 else if((string)inp.GetValue("Discontinuity Capturing") == "DC-minimum") solpar.iDC = 3; 590 else { 591 cout<< "Condition not defined for Discontinuity Capturing \n "; 592 exit(1); 593 } 594 595 // SCALAR DISCONTINUITY CAPTURING 596 597 vector<int> ivec = inp.GetValue("Scalar Discontinuity Capturing"); 598 for(i=0; i< 2; i++) solpar.idcsclr[i] = ivec[i]; 599 ivec.erase(ivec.begin(),ivec.end()); 600 601 602 // if((string)inp.GetValue("Scalar Discontinuity Capturing") == "No") solpar.idcsclr = 0; 603 // else if((string)inp.GetValue("Scalar Discontinuity Capturing") == "1") solpar.idcsclr = 1; 604 // else if((string)inp.GetValue("Scalar Discontinuity Capturing") == "2") solpar.idcsclr = 2; 605 // else { 606 // cout<< "Condition not defined for Scalar Discontinuity Capturing \n "; 607 // exit(1); 608 // } 609 if((string)inp.GetValue("Include Viscous Correction in Stabilization") == "True") 610 { 611 if(genpar.ipord == 1 ) genpar.idiff = 1; 612 else genpar.idiff = 2; 613 } 614 else { genpar.idiff = 0;} 615 616 genpar.irampViscOutlet = (int)inp.GetValue("Ramp Up Viscosity Near Outlet"); 617 618 genpar.istretchOutlet = (int)inp.GetValue("Stretch X Coordinate Near Outlet"); 619 620 genpar.iremoveStabTimeTerm = (int)inp.GetValue("Remove Time Term from Stabilization"); 621 622 timdat.flmpl = inp.GetValue("Lumped Mass Fraction on Left-hand-side"); 623 timdat.flmpr = inp.GetValue("Lumped Mass Fraction on Right-hand-side"); 624 625 timdat.iCFLworst = 0; 626 if((string)inp.GetValue("Dump CFL") == "True") 627 timdat.iCFLworst = 1; 628 629 intdat.intg[0][0]=inp.GetValue("Quadrature Rule on Interior"); 630 intdat.intg[0][1]=inp.GetValue("Quadrature Rule on Boundary"); 631 genpar.ibksiz = inp.GetValue("Number of Elements Per Block"); 632 633 ((string)inp.GetValue("Turn Off Source Terms for Scalars") 634 == "True")? sclrs.nosource=1:sclrs.nosource=0; 635 sclrs.tdecay=inp.GetValue("Decay Multiplier for Scalars"); 636 637 // TURBULENCE MODELING PARAMETER 638 int tpturb = turbvari.iles-turbvari.irans; 639 int ifrule; 640 if( tpturb != 0 ){ 641 642 643 turbvari.nohomog =inp.GetValue("Number of Homogenous Directions"); 644 645 if((string)inp.GetValue("Turbulence Wall Model Type") == "Slip Velocity") turbvar.itwmod = 1; 646 else if((string)inp.GetValue("Turbulence Wall Model Type") == "Effective Viscosity") turbvar.itwmod = 2; 647 else turbvar.itwmod = 0; 648 if (turbvari.irans < 0) turbvar.itwmod = turbvar.itwmod*(-1); 649 ifrule = inp.GetValue("Velocity Averaging Steps"); 650 turbvar.wtavei =(ifrule >0)? 1.0/ifrule : -1.0/ifrule; 651 652 if(turbvari.iles == 1){ 653 654 if((string)inp.GetValue("Dynamic Model Type") == "Bardina") turbvari.iles += 10; 655 else if((string)inp.GetValue("Dynamic Model Type") == "Projection") turbvari.iles += 20; 656 657 ifrule = inp.GetValue("Filter Integration Rule"); 658 turbvari.iles += ifrule-1; 659 ifrule = inp.GetValue("Dynamic Model Averaging Steps"); 660 turbvar.dtavei = (ifrule >0)? 1.0/ifrule : -1.0/ifrule; 661 turbvar.fwr1 = inp.GetValue("Filter Width Ratio"); 662 turbvar.flump = inp.GetValue("Lumping Factor for Filter"); 663 664 665 if ((string)inp.GetValue("Model Statistics") == "True" ) { 666 turbvari.modlstats = 1; } 667 else { 668 turbvari.modlstats = 0; } 669 670 if ((string)inp.GetValue("Double Filter") == "True" ) { 671 turbvari.i2filt = 1; } 672 else { 673 turbvari.i2filt = 0; } 674 675 if ((string)inp.GetValue("Model/SUPG Dissipation") == "True" ) { 676 turbvari.idis = 1; } 677 else { 678 turbvari.idis = 0; } 679 680 681 if((string)inp.GetValue("Dynamic Model Type") == "Standard") { 682 683 if((string)inp.GetValue("Dynamic Sub-Model Type") == "None") 684 turbvari.isubmod = 0; 685 else if((string)inp.GetValue("Dynamic Sub-Model Type") =="DFWR") 686 turbvari.isubmod = 1; 687 else if((string)inp.GetValue("Dynamic Sub-Model Type") =="SUPG") 688 turbvari.isubmod = 2; 689 } 690 else if((string)inp.GetValue("Dynamic Model Type") == "Projection") { 691 692 if((string)inp.GetValue("Projection Filter Type") == "Linear") 693 turbvari.ifproj = 0; 694 else if((string)inp.GetValue("Projection Filter Type") =="Quadratic") 695 turbvari.ifproj = 1; 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") =="ConsistentProj") 700 turbvari.isubmod = 1; 701 } 702 703 } 704 } 705 706 // SPEBC MODELING PARAMETERS 707 708 if ( (spebcvr.irscale = inp.GetValue("SPEBC Model Active")) >= 0 ){ 709 710 ifrule = inp.GetValue("Velocity Averaging Steps"); 711 turbvar.wtavei =(ifrule >0)? 1.0/ifrule : 1.0/inpdat.nstep[0]; 712 spebcvr.intpres = inp.GetValue("Interpolate Pressure"); 713 spebcvr.plandist = inp.GetValue("Distance between Planes"); 714 spebcvr.thetag = inp.GetValue("Theta Angle of Arc"); 715 spebcvr.ds = inp.GetValue("Distance for Velocity Averaging"); 716 spebcvr.tolerence = inp.GetValue("SPEBC Cylindrical Tolerance"); 717 spebcvr.radcyl = inp.GetValue("Radius of recycle plane"); 718 spebcvr.rbltin = inp.GetValue("Inlet Boundary Layer Thickness"); 719 spebcvr.rvscal = inp.GetValue("Vertical Velocity Scale Factor"); 720 } 721 722 // CARDIOVASCULAR MODELING PARAMETERS 723 if ( (string)inp.GetValue("Time Varying Boundary Conditions From File") == "True") 724 nomodule.itvn = 1; 725 else 726 nomodule.itvn = 0; 727 728 if ( nomodule.itvn ==1) 729 nomodule.bcttimescale = inp.GetValue("BCT Time Scale Factor"); 730 731 nomodule.ipvsq=0; 732 if(nomodule.icardio = inp.GetValue("Number of Coupled Surfaces")){ 733 if ( nomodule.icardio > MAXSURF ) { 734 cout << "Number of Coupled Surfaces > MAXSURF \n"; 735 exit(1); 736 } 737 if ( (string)inp.GetValue("Pressure Coupling") == "None") 738 nomodule.ipvsq=0; 739 if ( (string)inp.GetValue("Pressure Coupling") == "Explicit") 740 nomodule.ipvsq=1; 741 if ( (string)inp.GetValue("Pressure Coupling") == "Implicit") 742 nomodule.ipvsq=2; 743 if ( (string)inp.GetValue("Pressure Coupling") == "P-Implicit") 744 nomodule.ipvsq=3; 745 746 if(nomodule.numResistSrfs=inp.GetValue("Number of Resistance Surfaces")){ 747 ivec = inp.GetValue("List of Resistance Surfaces"); 748 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistResist[i] = 0; 749 for(i=0; i< nomodule.numResistSrfs; i++){ 750 nomodule.nsrflistResist[i+1] = ivec[i]; 751 } 752 vec = inp.GetValue("Resistance Values"); 753 for(i =0; i< MAXSURF+1 ; i++) nomodule.ValueListResist[i] = 0; 754 for(i =0; i< nomodule.numResistSrfs ; i++) nomodule.ValueListResist[i+1] = vec[i]; 755 vec.erase(vec.begin(),vec.end()); 756 } 757 if(nomodule.numImpSrfs=inp.GetValue("Number of Impedance Surfaces")){ 758 ivec = inp.GetValue("List of Impedance Surfaces"); 759 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistImp[i] = 0; 760 for(i=0; i< nomodule.numImpSrfs; i++){ 761 nomodule.nsrflistImp[i+1] = ivec[i]; 762 } 763 if ( (string)inp.GetValue("Impedance From File") == "True") 764 nomodule.impfile = 1; else nomodule.impfile = 0; 765 } 766 if(nomodule.numRCRSrfs=inp.GetValue("Number of RCR Surfaces")){ 767 ivec = inp.GetValue("List of RCR Surfaces"); 768 for(i=0;i<MAXSURF+1; i++) nomodule.nsrflistRCR[i] = 0; 769 for(i=0; i< nomodule.numRCRSrfs; i++){ 770 nomodule.nsrflistRCR[i+1] = ivec[i]; 771 } 772 if ( (string)inp.GetValue("RCR Values From File") == "True") 773 nomodule.ircrfile = 1; else nomodule.ircrfile = 0; 774 } 775 } 776 nomodule.ideformwall = 0; 777 if((string)inp.GetValue("Deformable Wall")=="True"){ 778 nomodule.ideformwall = 1; 779 nomodule.rhovw = inp.GetValue("Density of Vessel Wall"); 780 nomodule.thicknessvw = inp.GetValue("Thickness of Vessel Wall"); 781 nomodule.evw = inp.GetValue("Young Mod of Vessel Wall"); 782 nomodule.rnuvw = inp.GetValue("Poisson Ratio of Vessel Wall"); 783 nomodule.rshearconstantvw = inp.GetValue("Shear Constant of Vessel Wall"); 784 if((string)inp.GetValue("Wall Mass Matrix for LHS") == "True") nomodule.iwallmassfactor = 1; 785 else nomodule.iwallmassfactor = 0; 786 if((string)inp.GetValue("Wall Stiffness Matrix for LHS") == "True") nomodule.iwallstiffactor = 1; 787 else nomodule.iwallstiffactor = 0; 788 } 789 nomodule.iviscflux = 1; 790 if((string)inp.GetValue("Viscous Flux Flag") == "True") nomodule.iviscflux = 1; 791 if((string)inp.GetValue("Viscous Flux Flag") == "False") nomodule.iviscflux = 0; 792 793 794 // Scaling Parameters Keywords 795 796 outpar.ro = inp.GetValue("Density"); 797 outpar.vel = inp.GetValue("Velocity"); 798 outpar.press = inp.GetValue("Pressure"); 799 outpar.temper = inp.GetValue("Temperature"); 800 outpar.entrop = inp.GetValue("Entropy"); 801 802 // Step Sequencing 803 804 805 ivec = inp.GetValue("Step Construction"); 806 sequence.seqsize = ivec.size(); 807 if( sequence.seqsize > 100 || sequence.seqsize < 2 ) 808 cerr<<"Sequence size must be between 2 and 100 "<<endl; 809 810 for(i=0; i< sequence.seqsize; i++) 811 sequence.stepseq[i] = ivec[i]; 812 813 } 814 catch ( exception &e ) { 815 cout << endl << "Input exception: " << e.what() << endl << endl; 816 ierr = 001; 817 print_error_code(ierr); 818 return ierr; 819 } 820 821 return ierr; 822 823 } 824 825 void print_error_code(int ierr) { 826 /* 827 Return Error codes: 828 0xx Input error 829 1xx Solution Control 830 105 Turbulence Model not supported 831 832 2xx Material Properties 833 834 3xx Output Control 835 836 4xx Discretization Control 837 838 5xx Scaling Parameters 839 840 6xx Linear Solver Control 841 601 linear solver type not supported 842 */ 843 cout << endl << endl << "Input error detected: " << endl << endl; 844 if ( ierr == 001 ) { 845 cout << endl << "Input Directive not understood" << endl << endl; 846 } 847 if ( ierr == 105 ) { 848 cout << endl << "Turbulence Model Not Supported" << endl << endl; 849 } 850 if ( ierr == 601 ) { 851 cout << endl << "Linear Solver Type Not Supported" << endl << endl; 852 } 853 854 } 855