1 static const char help[] = "Integrate chemistry using TChem.\n"; 2 3 #include <petscts.h> 4 5 #if defined(PETSC_HAVE_TCHEM) 6 #if defined(MAX) 7 #undef MAX 8 #endif 9 #if defined(MIN) 10 #undef MIN 11 #endif 12 # include <TC_params.h> 13 # include <TC_interface.h> 14 #else 15 # error TChem is required for this example. Reconfigure PETSc using --download-tchem. 16 #endif 17 /* 18 See extchem.example.1 for how to run an example 19 20 See also h2_10sp.inp for another example 21 22 Determine sensitivity of final tempature on each variables initial conditions 23 -ts_dt 1.e-5 -ts_type cn -ts_adjoint_solve -ts_adjoint_view_solution draw 24 25 The solution for component i = 0 is the temperature. 26 27 The solution, i > 0, is the mass fraction, massf[i], of species i, i.e. mass of species i/ total mass of all species 28 29 The mole fraction molef[i], i > 0, is the number of moles of a species/ total number of moles of all species 30 Define M[i] = mass per mole of species i then 31 molef[i] = massf[i]/(M[i]*(sum_j massf[j]/M[j])) 32 33 FormMoleFraction(User,massf,molef) converts the mass fraction solution of each species to the mole fraction of each species. 34 35 36 These are other data sets for other possible runs 37 https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/n_heptane_v3.1_therm.dat 38 https://www-pls.llnl.gov/data/docs/science_and_technology/chemistry/combustion/nc7_ver3.1_mech.txt 39 40 41 */ 42 typedef struct _User *User; 43 struct _User { 44 PetscReal pressure; 45 int Nspec; 46 int Nreac; 47 PetscReal Tini; 48 double *tchemwork; 49 double *Jdense; /* Dense array workspace where Tchem computes the Jacobian */ 50 PetscInt *rows; 51 char **snames; 52 }; 53 54 55 static PetscErrorCode PrintSpecies(User,Vec); 56 static PetscErrorCode MassFractionToMoleFraction(User,Vec,Vec*); 57 static PetscErrorCode MoleFractionToMassFraction(User,Vec,Vec*); 58 static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*); 59 static PetscErrorCode FormRHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 60 static PetscErrorCode FormInitialSolution(TS,Vec,void*); 61 static PetscErrorCode ComputeMassConservation(Vec,PetscReal*,void*); 62 static PetscErrorCode MonitorMassConservation(TS,PetscInt,PetscReal,Vec,void*); 63 static PetscErrorCode MonitorTempature(TS,PetscInt,PetscReal,Vec,void*); 64 65 #define TCCHKERRQ(ierr) do {if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in TChem library, return code %d",ierr);} while (0) 66 67 int main(int argc,char **argv) 68 { 69 TS ts; /* time integrator */ 70 TSAdapt adapt; 71 Vec X,lambda; /* solution vector */ 72 Mat J; /* Jacobian matrix */ 73 PetscInt steps; 74 PetscErrorCode ierr; 75 PetscReal ftime,dt; 76 char chemfile[PETSC_MAX_PATH_LEN],thermofile[PETSC_MAX_PATH_LEN],lchemfile[PETSC_MAX_PATH_LEN],lthermofile[PETSC_MAX_PATH_LEN],lperiodic[PETSC_MAX_PATH_LEN]; 77 const char *periodic = "file://${PETSC_DIR}/${PETSC_ARCH}/share/periodictable.dat"; 78 struct _User user; /* user-defined work context */ 79 TSConvergedReason reason; 80 char **snames,*names; 81 PetscInt i; 82 TSTrajectory tj; 83 PetscBool flg = PETSC_FALSE,tflg = PETSC_FALSE,found; 84 85 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 86 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Chemistry solver options","");CHKERRQ(ierr); 87 ierr = PetscOptionsString("-chem","CHEMKIN input file","",chemfile,chemfile,sizeof(chemfile),NULL);CHKERRQ(ierr); 88 ierr = PetscFileRetrieve(PETSC_COMM_WORLD,chemfile,lchemfile,PETSC_MAX_PATH_LEN,&found);CHKERRQ(ierr); 89 if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot download %s and no local version %s",chemfile,lchemfile); 90 ierr = PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL);CHKERRQ(ierr); 91 ierr = PetscFileRetrieve(PETSC_COMM_WORLD,thermofile,lthermofile,PETSC_MAX_PATH_LEN,&found);CHKERRQ(ierr); 92 if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot download %s and no local version %s",thermofile,lthermofile); 93 user.pressure = 1.01325e5; /* Pascal */ 94 ierr = PetscOptionsReal("-pressure","Pressure of reaction [Pa]","",user.pressure,&user.pressure,NULL);CHKERRQ(ierr); 95 user.Tini = 1000; /* Kelvin */ 96 ierr = PetscOptionsReal("-Tini","Initial temperature [K]","",user.Tini,&user.Tini,NULL);CHKERRQ(ierr); 97 ierr = PetscOptionsBool("-monitor_mass","Monitor the total mass at each timestep","",flg,&flg,NULL);CHKERRQ(ierr); 98 ierr = PetscOptionsBool("-monitor_temp","Monitor the tempature each timestep","",tflg,&tflg,NULL);CHKERRQ(ierr); 99 ierr = PetscOptionsEnd();CHKERRQ(ierr); 100 101 /* tchem requires periodic table in current directory */ 102 ierr = PetscFileRetrieve(PETSC_COMM_WORLD,periodic,lperiodic,PETSC_MAX_PATH_LEN,&found);CHKERRQ(ierr); 103 if (!found) SETERRQ2(PETSC_COMM_WORLD,PETSC_ERR_FILE_OPEN,"Cannot located required periodic table %s or local version %s",periodic,lperiodic); 104 105 ierr = TC_initChem(lchemfile, lthermofile, 0, 1.0);TCCHKERRQ(ierr); 106 TC_setThermoPres(user.pressure); 107 user.Nspec = TC_getNspec(); 108 user.Nreac = TC_getNreac(); 109 /* 110 Get names of all species in easy to use array 111 */ 112 ierr = PetscMalloc1((user.Nspec+1)*LENGTHOFSPECNAME,&names);CHKERRQ(ierr); 113 ierr = PetscStrcpy(names,"Temp");CHKERRQ(ierr); 114 TC_getSnames(user.Nspec,names+LENGTHOFSPECNAME);CHKERRQ(ierr); 115 ierr = PetscMalloc1((user.Nspec+2),&snames);CHKERRQ(ierr); 116 for (i=0; i<user.Nspec+1; i++) snames[i] = names+i*LENGTHOFSPECNAME; 117 snames[user.Nspec+1] = NULL; 118 ierr = PetscStrArrayallocpy((const char *const *)snames,&user.snames);CHKERRQ(ierr); 119 ierr = PetscFree(snames);CHKERRQ(ierr); 120 ierr = PetscFree(names);CHKERRQ(ierr); 121 122 ierr = PetscMalloc3(user.Nspec+1,&user.tchemwork,PetscSqr(user.Nspec+1),&user.Jdense,user.Nspec+1,&user.rows);CHKERRQ(ierr); 123 ierr = VecCreateSeq(PETSC_COMM_SELF,user.Nspec+1,&X);CHKERRQ(ierr); 124 125 /* ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,PETSC_DECIDE,NULL,&J);CHKERRQ(ierr); */ 126 ierr = MatCreateSeqDense(PETSC_COMM_SELF,user.Nspec+1,user.Nspec+1,NULL,&J);CHKERRQ(ierr); 127 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 128 129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130 Create timestepping solver context 131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 132 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 133 ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr); 134 ierr = TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);CHKERRQ(ierr); 135 ierr = TSARKIMEXSetType(ts,TSARKIMEX4);CHKERRQ(ierr); 136 ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr); 137 ierr = TSSetRHSJacobian(ts,J,J,FormRHSJacobian,&user);CHKERRQ(ierr); 138 139 if (flg) { 140 ierr = TSMonitorSet(ts,MonitorMassConservation,NULL,NULL);CHKERRQ(ierr); 141 } 142 if (tflg) { 143 ierr = TSMonitorSet(ts,MonitorTempature,&user,NULL);CHKERRQ(ierr); 144 } 145 146 ftime = 1.0; 147 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr); 148 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 149 150 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 151 Set initial conditions 152 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 153 ierr = FormInitialSolution(ts,X,&user);CHKERRQ(ierr); 154 ierr = TSSetSolution(ts,X);CHKERRQ(ierr); 155 dt = 1e-10; /* Initial time step */ 156 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 157 ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr); 158 ierr = TSAdaptSetStepLimits(adapt,1e-12,1e-4);CHKERRQ(ierr); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */ 159 ierr = TSSetMaxSNESFailures(ts,-1);CHKERRQ(ierr); /* Retry step an unlimited number of times */ 160 161 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 162 Set runtime options 163 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 164 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 165 166 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 167 Set final conditions for sensitivities 168 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 169 ierr = VecDuplicate(X,&lambda);CHKERRQ(ierr); 170 ierr = TSSetCostGradients(ts,1,&lambda,NULL);CHKERRQ(ierr); 171 ierr = VecSetValue(lambda,0,1.0,INSERT_VALUES);CHKERRQ(ierr); 172 ierr = VecAssemblyBegin(lambda);CHKERRQ(ierr); 173 ierr = VecAssemblyEnd(lambda);CHKERRQ(ierr); 174 175 ierr = TSGetTrajectory(ts,&tj);CHKERRQ(ierr); 176 if (tj) { 177 ierr = TSTrajectorySetVariableNames(tj,(const char * const *)user.snames);CHKERRQ(ierr); 178 ierr = TSTrajectorySetTransform(tj,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user);CHKERRQ(ierr); 179 } 180 181 182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 183 Pass information to graphical monitoring routine 184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 185 ierr = TSMonitorLGSetVariableNames(ts,(const char * const *)user.snames);CHKERRQ(ierr); 186 ierr = TSMonitorLGSetTransform(ts,(PetscErrorCode (*)(void*,Vec,Vec*))MassFractionToMoleFraction,NULL,&user);CHKERRQ(ierr); 187 188 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 189 Solve ODE 190 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 191 ierr = TSSolve(ts,X);CHKERRQ(ierr); 192 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 193 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 194 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 195 ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);CHKERRQ(ierr); 196 197 /* { 198 Vec max; 199 PetscInt i; 200 const PetscReal *bmax; 201 202 ierr = TSMonitorEnvelopeGetBounds(ts,&max,NULL);CHKERRQ(ierr); 203 if (max) { 204 ierr = VecGetArrayRead(max,&bmax);CHKERRQ(ierr); 205 ierr = PetscPrintf(PETSC_COMM_SELF,"Species - maximum mass fraction\n");CHKERRQ(ierr); 206 for (i=1; i<user.Nspec; i++) { 207 if (bmax[i] > .01) {ierr = PetscPrintf(PETSC_COMM_SELF,"%s %g\n",user.snames[i],(double)bmax[i]);CHKERRQ(ierr);} 208 } 209 ierr = VecRestoreArrayRead(max,&bmax);CHKERRQ(ierr); 210 } 211 } 212 213 Vec y; 214 ierr = MassFractionToMoleFraction(&user,X,&y);CHKERRQ(ierr); 215 ierr = PrintSpecies(&user,y);CHKERRQ(ierr); 216 ierr = VecDestroy(&y);CHKERRQ(ierr); */ 217 218 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 219 Free work space. 220 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 221 TC_reset(); 222 ierr = PetscStrArrayDestroy(&user.snames);CHKERRQ(ierr); 223 ierr = MatDestroy(&J);CHKERRQ(ierr); 224 ierr = VecDestroy(&X);CHKERRQ(ierr); 225 ierr = VecDestroy(&lambda);CHKERRQ(ierr); 226 ierr = TSDestroy(&ts);CHKERRQ(ierr); 227 ierr = PetscFree3(user.tchemwork,user.Jdense,user.rows);CHKERRQ(ierr); 228 ierr = PetscFinalize(); 229 return ierr; 230 } 231 232 static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr) 233 { 234 User user = (User)ptr; 235 PetscErrorCode ierr; 236 PetscScalar *f; 237 const PetscScalar *x; 238 239 PetscFunctionBeginUser; 240 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 241 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 242 243 ierr = PetscArraycpy(user->tchemwork,x,user->Nspec+1);CHKERRQ(ierr); 244 user->tchemwork[0] *= user->Tini; /* Dimensionalize */ 245 ierr = TC_getSrc(user->tchemwork,user->Nspec+1,f);TCCHKERRQ(ierr); 246 f[0] /= user->Tini; /* Non-dimensionalize */ 247 248 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 249 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 static PetscErrorCode FormRHSJacobian(TS ts,PetscReal t,Vec X,Mat Amat,Mat Pmat,void *ptr) 254 { 255 User user = (User)ptr; 256 PetscErrorCode ierr; 257 const PetscScalar *x; 258 PetscInt M = user->Nspec+1,i; 259 260 PetscFunctionBeginUser; 261 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 262 ierr = PetscArraycpy(user->tchemwork,x,user->Nspec+1);CHKERRQ(ierr); 263 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 264 user->tchemwork[0] *= user->Tini; /* Dimensionalize temperature (first row) because that is what Tchem wants */ 265 ierr = TC_getJacTYN(user->tchemwork,user->Nspec,user->Jdense,1);CHKERRQ(ierr); 266 267 for (i=0; i<M; i++) user->Jdense[i + 0*M] /= user->Tini; /* Non-dimensionalize first column */ 268 for (i=0; i<M; i++) user->Jdense[0 + i*M] /= user->Tini; /* Non-dimensionalize first row */ 269 for (i=0; i<M; i++) user->rows[i] = i; 270 ierr = MatSetOption(Pmat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 271 ierr = MatSetOption(Pmat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 272 ierr = MatZeroEntries(Pmat);CHKERRQ(ierr); 273 ierr = MatSetValues(Pmat,M,user->rows,M,user->rows,user->Jdense,INSERT_VALUES);CHKERRQ(ierr); 274 275 ierr = MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 276 ierr = MatAssemblyEnd(Pmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 277 if (Amat != Pmat) { 278 ierr = MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 279 ierr = MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 280 } 281 PetscFunctionReturn(0); 282 } 283 284 PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx) 285 { 286 PetscScalar *x; 287 PetscErrorCode ierr; 288 PetscInt i; 289 Vec y; 290 const PetscInt maxspecies = 10; 291 PetscInt smax = maxspecies,mmax = maxspecies; 292 char *names[maxspecies]; 293 PetscReal molefracs[maxspecies],sum; 294 PetscBool flg; 295 296 PetscFunctionBeginUser; 297 ierr = VecZeroEntries(X);CHKERRQ(ierr); 298 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 299 x[0] = 1.0; /* Non-dimensionalized by user->Tini */ 300 301 ierr = PetscOptionsGetStringArray(NULL,NULL,"-initial_species",names,&smax,&flg);CHKERRQ(ierr); 302 if (smax < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Must provide at least two initial species"); 303 ierr = PetscOptionsGetRealArray(NULL,NULL,"-initial_mole",molefracs,&mmax,&flg);CHKERRQ(ierr); 304 if (smax != mmax) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Must provide same number of initial species %D as initial moles %D",smax,mmax); 305 sum = 0; 306 for (i=0; i<smax; i++) sum += molefracs[i]; 307 for (i=0; i<smax; i++) molefracs[i] = molefracs[i]/sum; 308 for (i=0; i<smax; i++) { 309 int ispec = TC_getSpos(names[i], strlen(names[i])); 310 if (ispec < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Could not find species %s",names[i]); 311 ierr = PetscPrintf(PETSC_COMM_SELF,"Species %d: %s %g\n",i,names[i],molefracs[i]);CHKERRQ(ierr); 312 x[1+ispec] = molefracs[i]; 313 } 314 for (i=0; i<smax; i++) { 315 ierr = PetscFree(names[i]);CHKERRQ(ierr); 316 } 317 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 318 /* PrintSpecies((User)ctx,X);CHKERRQ(ierr); */ 319 ierr = MoleFractionToMassFraction((User)ctx,X,&y);CHKERRQ(ierr); 320 ierr = VecCopy(y,X);CHKERRQ(ierr); 321 ierr = VecDestroy(&y);CHKERRQ(ierr); 322 PetscFunctionReturn(0); 323 } 324 325 /* 326 Converts the input vector which is in mass fractions (used by tchem) to mole fractions 327 */ 328 PetscErrorCode MassFractionToMoleFraction(User user,Vec massf,Vec *molef) 329 { 330 PetscErrorCode ierr; 331 PetscScalar *mof; 332 const PetscScalar *maf; 333 334 PetscFunctionBegin; 335 ierr = VecDuplicate(massf,molef);CHKERRQ(ierr); 336 ierr = VecGetArrayRead(massf,&maf);CHKERRQ(ierr); 337 ierr = VecGetArray(*molef,&mof);CHKERRQ(ierr); 338 mof[0] = maf[0]; /* copy over temperature */ 339 TC_getMs2Ml((double*)maf+1,user->Nspec,mof+1); 340 ierr = VecRestoreArray(*molef,&mof);CHKERRQ(ierr); 341 ierr = VecRestoreArrayRead(massf,&maf);CHKERRQ(ierr); 342 PetscFunctionReturn(0); 343 } 344 345 /* 346 Converts the input vector which is in mole fractions to mass fractions (used by tchem) 347 */ 348 PetscErrorCode MoleFractionToMassFraction(User user,Vec molef,Vec *massf) 349 { 350 PetscErrorCode ierr; 351 const PetscScalar *mof; 352 PetscScalar *maf; 353 354 PetscFunctionBegin; 355 ierr = VecDuplicate(molef,massf);CHKERRQ(ierr); 356 ierr = VecGetArrayRead(molef,&mof);CHKERRQ(ierr); 357 ierr = VecGetArray(*massf,&maf);CHKERRQ(ierr); 358 maf[0] = mof[0]; /* copy over temperature */ 359 TC_getMl2Ms((double*)mof+1,user->Nspec,maf+1); 360 ierr = VecRestoreArrayRead(molef,&mof);CHKERRQ(ierr); 361 ierr = VecRestoreArray(*massf,&maf);CHKERRQ(ierr); 362 PetscFunctionReturn(0); 363 } 364 365 PetscErrorCode ComputeMassConservation(Vec x,PetscReal *mass,void* ctx) 366 { 367 PetscErrorCode ierr; 368 369 PetscFunctionBegin; 370 ierr = VecSum(x,mass);CHKERRQ(ierr); 371 PetscFunctionReturn(0); 372 } 373 374 PetscErrorCode MonitorMassConservation(TS ts,PetscInt step,PetscReal time,Vec x,void* ctx) 375 { 376 const PetscScalar *T; 377 PetscReal mass; 378 PetscErrorCode ierr; 379 380 PetscFunctionBegin; 381 ierr = ComputeMassConservation(x,&mass,ctx);CHKERRQ(ierr); 382 ierr = VecGetArrayRead(x,&T);CHKERRQ(ierr); 383 mass -= PetscAbsScalar(T[0]); 384 ierr = VecRestoreArrayRead(x,&T);CHKERRQ(ierr); 385 ierr = PetscPrintf(PETSC_COMM_WORLD,"Timestep %D time %g percent mass lost or gained %g\n",step,(double)time,(double)100.*(1.0 - mass));CHKERRQ(ierr); 386 PetscFunctionReturn(0); 387 } 388 389 PetscErrorCode MonitorTempature(TS ts,PetscInt step,PetscReal time,Vec x,void* ctx) 390 { 391 User user = (User) ctx; 392 const PetscScalar *T; 393 PetscErrorCode ierr; 394 395 PetscFunctionBegin; 396 ierr = VecGetArrayRead(x,&T);CHKERRQ(ierr); 397 ierr = PetscPrintf(PETSC_COMM_WORLD,"Timestep %D time %g tempature %g\n",step,(double)time,(double)T[0]*user->Tini);CHKERRQ(ierr); 398 ierr = VecRestoreArrayRead(x,&T);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 402 /* 403 Prints out each species with its name 404 */ 405 PETSC_UNUSED PetscErrorCode PrintSpecies(User user,Vec molef) 406 { 407 PetscErrorCode ierr; 408 const PetscScalar *mof; 409 PetscInt i,*idx,n = user->Nspec+1; 410 411 PetscFunctionBegin; 412 ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr); 413 for (i=0; i<n;i++) idx[i] = i; 414 ierr = VecGetArrayRead(molef,&mof);CHKERRQ(ierr); 415 ierr = PetscSortRealWithPermutation(n,mof,idx);CHKERRQ(ierr); 416 for (i=0; i<n; i++) { 417 ierr = PetscPrintf(PETSC_COMM_SELF,"%6s %g\n",user->snames[idx[n-i-1]],mof[idx[n-i-1]]);CHKERRQ(ierr); 418 } 419 ierr = PetscFree(idx);CHKERRQ(ierr); 420 ierr = VecRestoreArrayRead(molef,&mof);CHKERRQ(ierr); 421 PetscFunctionReturn(0); 422 } 423 424 /*TEST 425 build: 426 requires: tchem 427 428 test: 429 args: -chem http://combustion.berkeley.edu/gri_mech/version30/files30/grimech30.dat -thermo http://combustion.berkeley.edu/gri_mech/version30/files30/thermo30.dat -initial_species CH4,O2,N2,AR -initial_mole 0.0948178320887,0.189635664177,0.706766236705,0.00878026702874 -Tini 1500 -Tini 1500 -ts_arkimex_fully_implicit -ts_max_snes_failures -1 -ts_adapt_dt_max 1e-4 -ts_arkimex_type 4 -ts_max_time .005 430 requires: !single 431 filter: grep -v iterations 432 433 TEST*/ 434