1 | /*************************************** 2 | $Header: /cvsroot/petscgraphics/chts.c,v 1.27 2004/03/07 01:56:18 hazelsct Exp $ 3 | 4 | This is the Cahn Hilliard timestepping code. It is provided here as an 5 | example usage of the Illuminator Distributed Visualization Library. 6 | ***************************************/ 7 | 8 | static char help[] = "Solves a nonlinear system in parallel using PETSc timestepping routines.\n\ 9 | \n\ 10 | The 2D or 3D transient Cahn-Hilliard phase field equation is solved on a 1x1\n\ 11 | square or 1x1x1 cube.\n\ 12 | The model is governed by the following parameters:\n\ 13 | -twodee : obvious (default is 3-D)\n\ 14 | -random : random initial condition (default is a box)\n\ 15 | -kappa <kap>, where <kap> = diffusivity\n\ 16 | -epsilon <eps>, where <eps> = interface thickness (default 1/mx)\n\ 17 | -gamma <gam>, where <gam> = interfacial tension (default 1)\n\ 18 | -mparam <m>, where <m> = free energy parameter, bet 0 & 1/2 for 0 stable\n\ 19 | Model parameters alpha and beta are set from epsilon and gamma according to:\n\ 20 | alpha = gam*eps beta=gam/eps\n\ 21 | Low-level options:\n\ 22 | -mx <xg>, where <xg> = number of grid points in the x-direction\n\ 23 | -my <yg>, where <yg> = number of grid points in the y-direction\n\ 24 | -mz <zg>, where <zg> = number of grid points in the z-direction\n\ 25 | -printg : print grid information\n\ 26 | Graphics of the contours of C are drawn by default at each timestep:\n\ 27 | -no_contours : do not draw contour plots of solution\n\ 28 | Parallelism can be invoked based on the DA construct:\n\ 29 | -Nx <npx>, where <npx> = number of processors in the x-direction\n\ 30 | -Ny <npy>, where <npy> = number of processors in the y-direction\n\ 31 | -Nz <npz>, where <npz> = number of processors in the z-direction\n\n"; 32 | 33 | 34 | /* Including cahnhill.h includes the necessary PETSc header files */ 35 | #include "cahnhill.h" 36 | #include "illuminator.h" 37 | 38 | 39 | /* Set #if to 1 to debug, 0 otherwise */ 40 | #undef DPRINTF 41 | #ifdef DEBUG 42 | #define DPRINTF(fmt, args...) ierr = PetscPrintf (PETSC_COMM_WORLD, "%s: " fmt, __FUNCT__, args); CHKERRQ(ierr) 43 | #else 44 | #define DPRINTF(fmt, args...) 45 | #endif 46 | 47 | 48 | /* Routines given below in this file */ 49 | int FormInitialCondition(AppCtx*,Vec); 50 | int UserLevelEnd(AppCtx*,Vec); 51 | int InitializeProblem(AppCtx*,Vec*); 52 | 53 | 54 | /*++++++++++++++++++++++++++++++++++++++ 55 | Wrapper for 56 | +latex+{\tt ch\_residual\_vector\_2d()} in {\tt cahnhill.c}. 57 | +html+ <tt>ch_residual_vector_2d()</tt> in <tt>cahnhill.c</tt>. 58 | 59 | int ch_ts_residual_vector_2d Usual return: zero or error. 60 | 61 | TS thets Timestepping context, ignored here. 62 | 63 | PetscScalar time Current time, also ignored. 64 | 65 | Vec X Current solution vector. 66 | 67 | Vec F Function vector to return. 68 | 69 | void *ptr User data pointer. 70 | ++++++++++++++++++++++++++++++++++++++*/ 71 | 72 | int ch_ts_residual_vector_2d 73 | (TS thets, PetscScalar time, Vec X, Vec F, void *ptr) 74 | { return ch_residual_vector_2d (X, F, ptr); } 75 | 76 | 77 | /*++++++++++++++++++++++++++++++++++++++ 78 | Wrapper for 79 | +latex+{\tt ch\_residual\_vector\_3d()} in {\tt cahnhill.c}. 80 | +html+ <tt>ch_residual_vector_3d()</tt> in <tt>cahnhill.c</tt>. 81 | 82 | int ch_ts_residual_vector_3d Usual return: zero or error. 83 | 84 | TS thets Timestepping context, ignored here. 85 | 86 | PetscScalar time Current time, also ignored. 87 | 88 | Vec X Current solution vector. 89 | 90 | Vec F Function vector to return. 91 | 92 | void *ptr User data pointer. 93 | ++++++++++++++++++++++++++++++++++++++*/ 94 | 95 | int ch_ts_residual_vector_3d 96 | (TS thets, PetscScalar time, Vec X, Vec F, void *ptr) 97 | { return ch_residual_vector_3d (X, F, ptr); } 98 | 99 | 100 | /*++++++++++++++++++++++++++++++++++++++ 101 | Wrapper for 102 | +latex+{\tt ch\_jacobian\_2d()} in {\tt cahnhill.c}. 103 | +html+ <tt>ch_jacobian_2d()</tt> in <tt>cahnhill.c</tt>. 104 | 105 | int ch_ts_jacobian_2d Usual return: zero or error. 106 | 107 | TS thets Timestepping context, ignored here. 108 | 109 | PetscScalar time Current time, also ignored. 110 | 111 | Vec X Current solution vector. 112 | 113 | Mat *A Place to put the new Jacobian. 114 | 115 | Mat *B Place to put the new conditioning matrix. 116 | 117 | MatStructure *flag Flag describing the volatility of the structure. 118 | 119 | void *ptr User data pointer. 120 | ++++++++++++++++++++++++++++++++++++++*/ 121 | 122 | int ch_ts_jacobian_2d (TS thets, PetscScalar time, Vec X, Mat *A, Mat *B, 123 | MatStructure *flag, void *ptr) { 124 | return ch_jacobian_2d (X, A, B, flag, ptr); } 125 | 126 | 127 | /*++++++++++++++++++++++++++++++++++++++ 128 | Wrapper for 129 | +latex+{\tt ch\_jacobian\_3d()} in {\tt cahnhill.c}. 130 | +html+ <tt>ch_jacobian_3d()</tt> in <tt>cahnhill.c</tt>. 131 | 132 | int ch_ts_jacobian_3d Usual return: zero or error. 133 | 134 | TS thets Timestepping context, ignored here. 135 | 136 | PetscScalar time Current time, also ignored. 137 | 138 | Vec X Current solution vector. 139 | 140 | Mat *A Place to put the new Jacobian. 141 | 142 | Mat *B Place to put the new conditioning matrix. 143 | 144 | MatStructure *flag Flag describing the volatility of the structure. 145 | 146 | void *ptr User data pointer. 147 | ++++++++++++++++++++++++++++++++++++++*/ 148 | 149 | int ch_ts_jacobian_3d (TS thets, PetscScalar time, Vec X, Mat *A, Mat *B, 150 | MatStructure *flag, void *ptr) { 151 | return ch_jacobian_3d (X, A, B, flag, ptr); } 152 | 153 | 154 | #undef __FUNCT__ 155 | #define __FUNCT__ "ch_ts_monitor" 156 | 157 | /*++++++++++++++++++++++++++++++++++++++ 158 | Monitor routine which displays the current state, using 159 | +latex+{\tt Illuminator}'s {\tt geomview} 160 | +html+ <tt>Illuminator</tt>'s <tt>geomview</tt> 161 | front-end (unless -no_contours is used); and also saves it using 162 | +latex+{\tt IlluMultiSave()} 163 | +html+ <tt>IlluMultiSave()</tt> 164 | if -save_data is specified. 165 | 166 | int ch_ts_monitor Usual return: zero or error. 167 | 168 | TS thets Timestepping context, ignored here. 169 | 170 | int stepno Current time step number. 171 | 172 | PetscScalar time Current time. 173 | 174 | Vec X Vector of current solved field values. 175 | 176 | void *ptr User data pointer. 177 | ++++++++++++++++++++++++++++++++++++++*/ 178 | 179 | int ch_ts_monitor (TS thets, int stepno, PetscScalar time, Vec X, void *ptr) 180 | { 181 | AppCtx *user; 182 | int temp, ierr; 183 | PetscReal xmin, xmax; 184 | PetscScalar minmax[6] = { 0.,1., 0.,1., 0.,1. }; 185 | /* Example colors and isoquant values to pass to DATriangulate() */ 186 | /* PetscScalar colors[16] = { 1,0,0,.5, 1,1,0,.5, 0,1,0,.5, 0,0,1,.5 }; 187 | PetscScalar isovals[4] = { .2, .4, .6, .8 }; */ 188 | 189 | user = (AppCtx *)ptr; 190 | 191 | ierr = VecMin (X, &temp, &xmin); CHKERRQ (ierr); 192 | ierr = VecMax (X, &temp, &xmax); CHKERRQ (ierr); 193 | PetscPrintf (user->comm,"Step %d, Time %g, C values from %g to %g\n", 194 | stepno, time, xmin, xmax); 195 | 196 | if (!user->no_contours) 197 | { 198 | if (user->threedee) 199 | { 200 | DPRINTF ("Starting triangulation\n",0); 201 | ierr = DATriangulate (user->da, X, user->chvar, minmax, PETSC_DECIDE, 202 | PETSC_NULL, PETSC_NULL, PETSC_FALSE, 203 | PETSC_FALSE, PETSC_FALSE); CHKERRQ (ierr); 204 | DPRINTF ("Done, sending to Geomview\n",0); 205 | ierr = GeomviewDisplayTriangulation 206 | (user->comm, minmax, "Cahn-Hilliard", PETSC_TRUE); CHKERRQ (ierr); 207 | DPRINTF ("Done.\n",0); 208 | } 209 | else 210 | { 211 | ierr = VecView (X,user->theviewer); CHKERRQ (ierr); 212 | } 213 | } 214 | 215 | if (user->save_data) 216 | { 217 | PetscReal deltat; 218 | field_plot_type chtype=FIELD_SCALAR_CONTOURS; 219 | char filename [100], *paramdata [4], 220 | *paramnames [4] = { "time", "timestep", "gamma", "kappa" }; 221 | 222 | for (ierr=0; ierr<4; ierr++) 223 | paramdata[ierr] = (char *) malloc (50*sizeof(char)); 224 | snprintf (filename, 99, "chtsout.time%.3d", stepno); 225 | TSGetTimeStep (thets, &deltat); 226 | snprintf (paramdata[0], 49, "%g", time); 227 | snprintf (paramdata[1], 49, "%g", deltat); 228 | snprintf (paramdata[2], 49, "%g", user->gamma); 229 | snprintf (paramdata[3], 49, "%g", user->kappa); 230 | DPRINTF ("Storing data using IlluMultiSave()\n",0); 231 | IlluMultiSave (user->da, X, filename, 1.,1.,1., &chtype, 4, paramnames, 232 | paramdata, COMPRESS_INT_NONE | COMPRESS_GZIP_FAST); 233 | } 234 | 235 | DPRINTF ("Completed timestep monitor tasks.\n",0); 236 | 237 | return 0; 238 | } 239 | 240 | 241 | #undef __FUNCT__ 242 | #define __FUNCT__ "main" 243 | 244 | /*++++++++++++++++++++++++++++++++++++++ 245 | The usual main function. 246 | 247 | int main Returns 0 or error. 248 | 249 | int argc Number of args. 250 | 251 | char **argv The args. 252 | ++++++++++++++++++++++++++++++++++++++*/ 253 | 254 | int main (int argc, char **argv) 255 | { 256 | TS thets; /* the timestepping solver */ 257 | Vec x; /* solution vector */ 258 | AppCtx user; /* user-defined work context */ 259 | int dim, ierr; 260 | PetscDraw draw; 261 | PetscTruth matfree; 262 | PetscReal xmin, xmax; 263 | int temp; 264 | PetscScalar ftime, time_total_max = 100.0; /* default max total time */ 265 | int steps = 0, time_steps_max = 5; /* default max timesteps */ 266 | 267 | PetscInitialize (&argc, &argv, (char *)0, help); 268 | ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &user.rank); CHKERRQ (ierr); 269 | ierr = MPI_Comm_size (PETSC_COMM_WORLD, &user.size); CHKERRQ (ierr); 270 | user.comm = PETSC_COMM_WORLD; 271 | user.mc = 1; user.chvar = 0; /* Sets number of variables, which is conc */ 272 | 273 | /* Create user context, set problem data, create vector data structures. 274 | Also, set the initial condition */ 275 | 276 | DPRINTF ("About to initialize problem\n",0); 277 | ierr = InitializeProblem (&user, &x); CHKERRQ (ierr); 278 | if (user.load_data > -1) 279 | steps = user.load_data; 280 | ierr = VecDuplicate (user.localX, &user.localF); CHKERRQ (ierr); 281 | 282 | /* Set up displays to show graph of the solution */ 283 | 284 | if (!user.no_contours) 285 | { 286 | if (user.threedee) 287 | { 288 | ierr = GeomviewBegin (PETSC_COMM_WORLD); CHKERRQ (ierr); 289 | } 290 | else 291 | { 292 | ierr = PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "", PETSC_DECIDE, 293 | PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 294 | &user.theviewer); CHKERRQ (ierr); 295 | ierr = PetscViewerDrawGetDraw (user.theviewer, 0, &draw); 296 | CHKERRQ (ierr); 297 | ierr = PetscDrawSetDoubleBuffer (draw); CHKERRQ (ierr); 298 | } 299 | } 300 | 301 | /* Create timestepping solver context */ 302 | DPRINTF ("Creating timestepping context\n",0); 303 | ierr = TSCreate (PETSC_COMM_WORLD, &thets); CHKERRQ (ierr); 304 | ierr = TSSetProblemType (thets, TS_NONLINEAR); CHKERRQ (ierr); 305 | ierr = VecGetSize (x, &dim); CHKERRQ (ierr); 306 | ierr = PetscPrintf (user.comm, "global size = %d, kappa = %g, epsilon = %g, " 307 | "gamma = %g, mparam = %g\n", dim, user.kappa, 308 | user.epsilon, user.gamma, user.mparam); CHKERRQ (ierr); 309 | 310 | /* Set function evaluation routine and monitor */ 311 | DPRINTF ("Setting RHS function\n",0); 312 | if (user.threedee) 313 | ierr = TSSetRHSFunction (thets, ch_ts_residual_vector_3d, &user); 314 | else 315 | ierr = TSSetRHSFunction (thets, ch_ts_residual_vector_2d, &user); 316 | CHKERRQ(ierr); 317 | ierr = TSSetMonitor (thets, ch_ts_monitor, &user, PETSC_NULL); CHKERRQ(ierr); 318 | 319 | /* This condition prevents the construction of the Jacobian if we're 320 | running matrix-free. */ 321 | ierr = PetscOptionsHasName (PETSC_NULL, "-snes_mf", &matfree); CHKERRQ(ierr); 322 | 323 | if (!matfree) 324 | { 325 | /* Set up the Jacobian using cahnhill.c subroutines */ 326 | DPRINTF ("Using analytical Jacobian from cahnhill.c\n",0); 327 | if (user.threedee) { 328 | ierr = ch_jacobian_alpha_3d (&user); CHKERRQ (ierr); 329 | ierr = TSSetRHSJacobian (thets, user.J, user.J, ch_ts_jacobian_3d, 330 | &user); CHKERRQ (ierr); } 331 | else { 332 | ierr = ch_jacobian_alpha_2d (&user); CHKERRQ (ierr); 333 | ierr = TSSetRHSJacobian (thets, user.J, user.J, ch_ts_jacobian_2d, 334 | &user); CHKERRQ (ierr);} 335 | /*}*/ 336 | } 337 | 338 | /* Set solution vector and initial timestep (currently a fraction of the 339 | explicit diffusion stability criterion */ 340 | ierr = TSSetInitialTimeStep (thets, 0.0, 0.1/(user.mx-1)/(user.mx-1)); 341 | CHKERRQ (ierr); 342 | ierr = TSSetSolution (thets, x); CHKERRQ (ierr); 343 | 344 | /* Customize timestepping solver, default to Crank-Nicholson */ 345 | ierr = TSSetDuration (thets, time_steps_max, time_total_max); CHKERRQ (ierr); 346 | ierr = TSSetType (thets, TS_CRANK_NICHOLSON); CHKERRQ (ierr); 347 | ierr = TSSetFromOptions (thets); CHKERRQ (ierr); 348 | 349 | /* Run the solver */ 350 | DPRINTF ("Running the solver\n",0); 351 | ierr = TSStep (thets, &steps, &ftime); CHKERRQ (ierr); 352 | 353 | /* Final clean-up */ 354 | DPRINTF ("Cleaning up\n",0); 355 | if (!user.no_contours) 356 | { 357 | if (user.threedee) 358 | { 359 | ierr = GeomviewEnd (PETSC_COMM_WORLD); CHKERRQ (ierr); 360 | } 361 | else 362 | { 363 | ierr = PetscViewerDestroy (user.theviewer); CHKERRQ (ierr); 364 | } 365 | } 366 | ierr = VecDestroy (user.localX); CHKERRQ (ierr); 367 | ierr = VecDestroy (x); CHKERRQ (ierr); 368 | ierr = VecDestroy (user.localF); CHKERRQ (ierr); 369 | ierr = TSDestroy (thets); CHKERRQ (ierr); 370 | ierr = PetscFree (user.label); CHKERRQ (ierr); 371 | ierr = DADestroy (user.da); CHKERRQ (ierr); 372 | 373 | PetscFinalize (); 374 | printf ("Game over man!\n"); 375 | return 0; 376 | } 377 | 378 | 379 | #undef __FUNCT__ 380 | #define __FUNCT__ "FormInitialCondition" 381 | 382 | /*++++++++++++++++++++++++++++++++++++++ 383 | Like it says, put together the initial condition. 384 | 385 | int FormInitialCondition Returns zero or error. 386 | 387 | AppCtx *user The user context structure. 388 | 389 | Vec X Vector in which to place the initial condition. 390 | ++++++++++++++++++++++++++++++++++++++*/ 391 | 392 | int FormInitialCondition (AppCtx *user, Vec X) 393 | { 394 | int i,j,k, row, mc, chvar, mx,my,mz, ierr, xs,ys,zs, xm,ym,zm; 395 | int gxm,gym,gzm, gxs,gys,gzs; 396 | PetscScalar mparam; 397 | PetscScalar *x; 398 | Vec localX = user->localX; 399 | PetscRandom rand; 400 | 401 | mc = user->mc; 402 | chvar = user->chvar; 403 | mx = user->mx; my = user->my; mz = user->mz; 404 | mparam = user->mparam; 405 | 406 | /* Get a pointer to vector data. 407 | - For default PETSc vectors, VecGetArray() returns a pointer to 408 | the data array. Otherwise, the routine is implementation dependent. 409 | - You MUST call VecRestoreArray() when you no longer need access to 410 | the array. */ 411 | if (user->random) 412 | { 413 | DPRINTF ("Setting initial condition as random numbers\n",0); 414 | ierr = PetscRandomCreate (user->comm, RANDOM_DEFAULT, &rand); 415 | CHKERRQ (ierr); 416 | ierr = PetscRandomSetInterval (rand, 0.35, 0.65); CHKERRQ (ierr); 417 | ierr = VecSetRandom (rand, X); CHKERRQ (ierr); 418 | ierr = PetscRandomDestroy (rand); CHKERRQ (ierr); 419 | } 420 | else if (user->load_data > -1) 421 | { 422 | int usermetacount; 423 | char basename [1000], **usermetanames, **usermetadata; 424 | sprintf (basename, "chtsout.time%.3d", user->load_data); 425 | DPRINTF ("Loading data for timestep %d, basename %s\n", user->load_data, 426 | basename); 427 | IlluMultiRead (user->da, X, basename, &usermetacount, &usermetanames, 428 | &usermetadata); 429 | /* TODO: free these */ 430 | for (i=0; i<usermetacount; i++) 431 | PetscPrintf (PETSC_COMM_WORLD, "Parameter \"%s\" = \"%s\"\n", 432 | usermetanames [i], usermetadata [i]); 433 | } 434 | else 435 | { 436 | DPRINTF ("Getting local array\n",0); 437 | ierr = VecGetArray(localX,&x); CHKERRQ (ierr); 438 | 439 | /* Get local grid boundaries (for 2-dimensional DA): 440 | xs, ys - starting grid indices (no ghost points) 441 | xm, ym - widths of local grid (no ghost points) 442 | gxs, gys - starting grid indices (including ghost points) 443 | gxm, gym - widths of local grid (including ghost points) */ 444 | DPRINTF ("Getting corners and ghost corners\n",0); 445 | ierr = DAGetCorners (user->da, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr); 446 | ierr = DAGetGhostCorners (user->da, &gxs,&gys,&gzs, &gxm,&gym,&gzm); 447 | CHKERRQ (ierr); 448 | 449 | /* Set up 2-D so it works */ 450 | if (!user->threedee) { zs = gzs = 0; zm = 1; mz = 10; } 451 | 452 | /* Compute initial condition over the locally owned part of the grid 453 | Initial condition is motionless fluid and equilibrium temperature */ 454 | DPRINTF ("Looping to set initial condition\n",0); 455 | for (k=zs; k<zs+zm; k++) 456 | for (j=ys; j<ys+ym; j++) 457 | for (i=xs; i<xs+xm; i++) 458 | { 459 | row = i - gxs + (j - gys)*gxm + (k-gzs)*gxm*gym; 460 | /* x[C(row)] = (i>=mx*0.43921) ? 1.0 : 0.0; */ 461 | x[C(row)] = ((i<mx/3 || i>2*mx/3) && (j<my/3 || j>2*my/3) && 462 | (k<mz/3 || k>2*mz/3)) ? 1.0 : 0.0; 463 | /* x[C(row)] = (i<mx/2 && j<my/2 && k<mz/2) ? 1.0 : 0.0; */ 464 | /* x[V(row)] = 0.0; 465 | x[Omega(row)] = 0.0; 466 | x[Temp(row)] = (mparam>0)*(PetscScalar)(i)/(PetscScalar)(mx-1); */ 467 | } 468 | 469 | /* Restore vector */ 470 | DPRINTF ("Restoring array to local vector\n",0); 471 | ierr = VecRestoreArray (localX,&x); CHKERRQ (ierr); 472 | 473 | /* Insert values into global vector */ 474 | DPRINTF ("Inserting local vector values into global vector\n",0); 475 | ierr = DALocalToGlobal (user->da,localX,INSERT_VALUES,X); CHKERRQ (ierr); 476 | } 477 | 478 | return 0; 479 | } 480 | 481 | 482 | #undef __FUNCT__ 483 | #define __FUNCT__ "InitializeProblem" 484 | 485 | /*++++++++++++++++++++++++++++++++++++++ 486 | This takes the gory details of initialization out of the way (importing 487 | parameters into the user context, etc.). 488 | 489 | int InitializeProblem Returns zero or error. 490 | 491 | AppCtx *user The user context to fill. 492 | 493 | Vec *xvec Vector into which to put the initial condition. 494 | ++++++++++++++++++++++++++++++++++++++*/ 495 | 496 | int InitializeProblem (AppCtx *user, Vec *xvec) 497 | { 498 | int Nx,Ny,Nz; /* number of processors in x-, y- and z- directions */ 499 | int xs,xm,ys,ym,zs,zm, Nlocal,ierr; 500 | Vec xv; 501 | PetscTruth twodee; 502 | 503 | /* Initialize problem parameters */ 504 | DPRINTF ("Initializing user->parameters\n",0); 505 | user->mx = 20; 506 | user->my = 20; 507 | user->mz = 20; 508 | ierr = PetscOptionsGetInt(PETSC_NULL, "-mx", &user->mx, PETSC_NULL); 509 | CHKERRQ (ierr); 510 | ierr = PetscOptionsGetInt(PETSC_NULL, "-my", &user->my, PETSC_NULL); 511 | CHKERRQ (ierr); 512 | ierr = PetscOptionsGetInt(PETSC_NULL, "-mz", &user->mz, PETSC_NULL); 513 | CHKERRQ (ierr); 514 | /* No. of components in the unknown vector and auxiliary vector */ 515 | user->mc = 1; 516 | /* Problem parameters (kappa, epsilon, gamma, and mparam) */ 517 | user->kappa = 1.0; 518 | ierr = PetscOptionsGetScalar(PETSC_NULL, "-kappa", &user->kappa, PETSC_NULL); 519 | CHKERRQ (ierr); 520 | user->epsilon = 1.0/(user->mx-1); 521 | ierr = PetscOptionsGetScalar (PETSC_NULL, "-epsilon", &user->epsilon, 522 | PETSC_NULL); CHKERRQ (ierr); 523 | user->gamma = 1.0; 524 | ierr = PetscOptionsGetScalar(PETSC_NULL, "-gamma", &user->gamma, PETSC_NULL); 525 | CHKERRQ (ierr); 526 | user->mparam = 0.0; 527 | ierr = PetscOptionsGetScalar (PETSC_NULL, "-mparam", &user->mparam, 528 | PETSC_NULL); CHKERRQ (ierr); 529 | ierr = PetscOptionsHasName (PETSC_NULL, "-twodee", &twodee); 530 | user->threedee = !twodee; 531 | CHKERRQ (ierr); 532 | ierr = PetscOptionsHasName (PETSC_NULL, "-printv", &user->print_vecs); 533 | CHKERRQ (ierr); 534 | ierr = PetscOptionsHasName (PETSC_NULL, "-printg", &user->print_grid); 535 | CHKERRQ (ierr); 536 | ierr = PetscOptionsHasName (PETSC_NULL, "-no_contours", &user->no_contours); 537 | CHKERRQ (ierr); 538 | ierr = PetscOptionsHasName (PETSC_NULL, "-random", &user->random); 539 | CHKERRQ (ierr); 540 | ierr = PetscOptionsHasName (PETSC_NULL, "-save_data", &user->save_data); 541 | CHKERRQ (ierr); 542 | user->load_data = -1; 543 | ierr = PetscOptionsGetInt (PETSC_NULL, "-load_data", &user->load_data, 544 | PETSC_NULL); CHKERRQ (ierr); 545 | 546 | /* Create distributed array (DA) to manage parallel grid and vectors 547 | for principal unknowns (x) and governing residuals (f) 548 | Note the stencil width of 2 for this 4th-order equation. */ 549 | DPRINTF ("Creating distributed array\n",0); 550 | Nx = PETSC_DECIDE; 551 | Ny = PETSC_DECIDE; 552 | Nz = PETSC_DECIDE; 553 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Nx", &Nx, PETSC_NULL);CHKERRQ(ierr); 554 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Ny", &Ny, PETSC_NULL);CHKERRQ(ierr); 555 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Nz", &Nz, PETSC_NULL);CHKERRQ(ierr); 556 | if (user->threedee) 557 | { 558 | DPRINTF ("Three Dee!\n",0); 559 | user->period = DA_XYZPERIODIC; 560 | ierr = DACreate3d (PETSC_COMM_WORLD, user->period, DA_STENCIL_BOX, 561 | user->mx, user->my, user->mz, Nx, Ny, Nz, user->mc, 2, 562 | PETSC_NULL, PETSC_NULL, PETSC_NULL, &user->da); 563 | CHKERRQ (ierr); 564 | } 565 | else 566 | { 567 | user->period = DA_XYPERIODIC; 568 | ierr = DACreate2d (PETSC_COMM_WORLD, user->period, DA_STENCIL_BOX, 569 | user->mx, user->my, Nx, Ny, user->mc, 2, 570 | PETSC_NULL, PETSC_NULL, &user->da); CHKERRQ (ierr); 571 | user->mz = Nz = 1; 572 | } 573 | 574 | /* Extract global vector from DA */ 575 | DPRINTF ("Extracting global vector from DA...\n",0); 576 | ierr = DACreateGlobalVector(user->da,&xv); CHKERRQ (ierr); 577 | 578 | /* Label PDE components */ 579 | DPRINTF ("Labeling PDE components\n",0); 580 | ierr = PetscMalloc (user->mc * sizeof(char*), &(user->label));CHKERRQ (ierr); 581 | user->label[0] = "Concentration (C)"; 582 | /* user->label[1] = "Velocity (V)"; 583 | user->label[2] = "Omega"; 584 | user->label[3] = "Temperature"; */ 585 | ierr = DASetFieldName (user->da, 0, user->label[0]); CHKERRQ(ierr); 586 | 587 | user->x_old = 0; 588 | 589 | /* Get local vector */ 590 | DPRINTF ("Getting local vector\n",0); 591 | ierr = DACreateLocalVector (user->da, &user->localX); CHKERRQ (ierr); 592 | 593 | /* Print grid info */ 594 | if (user->print_grid) 595 | { 596 | DPRINTF ("Printing grid information\n",0); 597 | ierr = DAView(user->da,PETSC_VIEWER_STDOUT_SELF); CHKERRQ (ierr); 598 | ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr); 599 | if (!user->threedee) { 600 | zs = 0; zm = 1; } 601 | ierr = PetscPrintf(PETSC_COMM_WORLD, 602 | "global grid: %d X %d X %d with %d components per" 603 | " node ==> global vector dimension %d\n", 604 | user->mx, user->my, user->mz, user->mc, 605 | user->mc*user->mx*user->my*user->mz); 606 | CHKERRQ (ierr); 607 | fflush(stdout); 608 | ierr = VecGetLocalSize (xv, &Nlocal); CHKERRQ (ierr); 609 | ierr = PetscSynchronizedPrintf 610 | (PETSC_COMM_WORLD,"[%d] local grid %d X %d X %d with %d components" 611 | " per node ==> local vector dimension %d\n", 612 | user->rank, xm, ym, zm, user->mc, Nlocal); CHKERRQ (ierr); 613 | ierr = PetscSynchronizedFlush (PETSC_COMM_WORLD); CHKERRQ (ierr); 614 | } 615 | 616 | /* Compute initial condition */ 617 | DPRINTF ("Forming inital condition\n",0); 618 | ierr = FormInitialCondition (user, xv); CHKERRQ (ierr); 619 | 620 | *xvec = xv; 621 | return 0; 622 | }