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