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