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