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