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