1 | /*************************************** 2 | $Header: /cvsroot/petscgraphics/chts.c,v 1.26 2003/07/25 15:41:51 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); CHKERRQ (ierr); 203 | DPRINTF ("Done, sending to Geomview\n",0); 204 | ierr = GeomviewDisplayTriangulation 205 | (user->comm, minmax, "Cahn-Hilliard", PETSC_TRUE); CHKERRQ (ierr); 206 | DPRINTF ("Done.\n",0); 207 | } 208 | else 209 | { 210 | ierr = VecView (X,user->theviewer); CHKERRQ (ierr); 211 | } 212 | } 213 | 214 | if (user->save_data) 215 | { 216 | PetscReal deltat; 217 | field_plot_type chtype=FIELD_SCALAR_CONTOURS; 218 | char filename [100], *paramdata [4], 219 | *paramnames [4] = { "time", "timestep", "gamma", "kappa" }; 220 | 221 | for (ierr=0; ierr<4; ierr++) 222 | paramdata[ierr] = (char *) malloc (50*sizeof(char)); 223 | snprintf (filename, 99, "chtsout.time%.3d", stepno); 224 | TSGetTimeStep (thets, &deltat); 225 | snprintf (paramdata[0], 49, "%g", time); 226 | snprintf (paramdata[1], 49, "%g", deltat); 227 | snprintf (paramdata[2], 49, "%g", user->gamma); 228 | snprintf (paramdata[3], 49, "%g", user->kappa); 229 | DPRINTF ("Storing data using IlluMultiSave()\n",0); 230 | IlluMultiSave (user->da, X, filename, 1.,1.,1., &chtype, 4, paramnames, 231 | paramdata, COMPRESS_INT_NONE | COMPRESS_GZIP_FAST); 232 | } 233 | 234 | DPRINTF ("Completed timestep monitor tasks.\n",0); 235 | 236 | return 0; 237 | } 238 | 239 | 240 | #undef __FUNCT__ 241 | #define __FUNCT__ "main" 242 | 243 | /*++++++++++++++++++++++++++++++++++++++ 244 | The usual main function. 245 | 246 | int main Returns 0 or error. 247 | 248 | int argc Number of args. 249 | 250 | char **argv The args. 251 | ++++++++++++++++++++++++++++++++++++++*/ 252 | 253 | int main (int argc, char **argv) 254 | { 255 | TS thets; /* the timestepping solver */ 256 | Vec x; /* solution vector */ 257 | AppCtx user; /* user-defined work context */ 258 | int dim, ierr; 259 | PetscDraw draw; 260 | PetscTruth matfree; 261 | PetscReal xmin, xmax; 262 | int temp; 263 | PetscScalar ftime, time_total_max = 100.0; /* default max total time */ 264 | int steps = 0, time_steps_max = 5; /* default max timesteps */ 265 | 266 | PetscInitialize (&argc, &argv, (char *)0, help); 267 | ierr = MPI_Comm_rank (PETSC_COMM_WORLD, &user.rank); CHKERRQ (ierr); 268 | ierr = MPI_Comm_size (PETSC_COMM_WORLD, &user.size); CHKERRQ (ierr); 269 | user.comm = PETSC_COMM_WORLD; 270 | user.mc = 1; user.chvar = 0; /* Sets number of variables, which is conc */ 271 | 272 | /* Create user context, set problem data, create vector data structures. 273 | Also, set the initial condition */ 274 | 275 | DPRINTF ("About to initialize problem\n",0); 276 | ierr = InitializeProblem (&user, &x); CHKERRQ (ierr); 277 | if (user.load_data > -1) 278 | steps = user.load_data; 279 | ierr = VecDuplicate (user.localX, &user.localF); CHKERRQ (ierr); 280 | 281 | /* Set up displays to show graph of the solution */ 282 | 283 | if (!user.no_contours) 284 | { 285 | if (user.threedee) 286 | { 287 | ierr = GeomviewBegin (PETSC_COMM_WORLD); CHKERRQ (ierr); 288 | } 289 | else 290 | { 291 | ierr = PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "", PETSC_DECIDE, 292 | PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 293 | &user.theviewer); CHKERRQ (ierr); 294 | ierr = PetscViewerDrawGetDraw (user.theviewer, 0, &draw); 295 | CHKERRQ (ierr); 296 | ierr = PetscDrawSetDoubleBuffer (draw); CHKERRQ (ierr); 297 | } 298 | } 299 | 300 | /* Create timestepping solver context */ 301 | DPRINTF ("Creating timestepping context\n",0); 302 | ierr = TSCreate (PETSC_COMM_WORLD, &thets); CHKERRQ (ierr); 303 | ierr = TSSetProblemType (thets, TS_NONLINEAR); CHKERRQ (ierr); 304 | ierr = VecGetSize (x, &dim); CHKERRQ (ierr); 305 | ierr = PetscPrintf (user.comm, "global size = %d, kappa = %g, epsilon = %g, " 306 | "gamma = %g, mparam = %g\n", dim, user.kappa, 307 | user.epsilon, user.gamma, user.mparam); CHKERRQ (ierr); 308 | 309 | /* Set function evaluation routine and monitor */ 310 | DPRINTF ("Setting RHS function\n",0); 311 | if (user.threedee) 312 | ierr = TSSetRHSFunction (thets, ch_ts_residual_vector_3d, &user); 313 | else 314 | ierr = TSSetRHSFunction (thets, ch_ts_residual_vector_2d, &user); 315 | CHKERRQ(ierr); 316 | ierr = TSSetMonitor (thets, ch_ts_monitor, &user, PETSC_NULL); CHKERRQ(ierr); 317 | 318 | /* This condition prevents the construction of the Jacobian if we're 319 | running matrix-free. */ 320 | ierr = PetscOptionsHasName (PETSC_NULL, "-snes_mf", &matfree); CHKERRQ(ierr); 321 | 322 | if (!matfree) 323 | { 324 | /* Set up the Jacobian using cahnhill.c subroutines */ 325 | DPRINTF ("Using analytical Jacobian from cahnhill.c\n",0); 326 | if (user.threedee) { 327 | ierr = ch_jacobian_alpha_3d (&user); CHKERRQ (ierr); 328 | ierr = TSSetRHSJacobian (thets, user.J, user.J, ch_ts_jacobian_3d, 329 | &user); CHKERRQ (ierr); } 330 | else { 331 | ierr = ch_jacobian_alpha_2d (&user); CHKERRQ (ierr); 332 | ierr = TSSetRHSJacobian (thets, user.J, user.J, ch_ts_jacobian_2d, 333 | &user); CHKERRQ (ierr);} 334 | /*}*/ 335 | } 336 | 337 | /* Set solution vector and initial timestep (currently a fraction of the 338 | explicit diffusion stability criterion */ 339 | ierr = TSSetInitialTimeStep (thets, 0.0, 0.1/(user.mx-1)/(user.mx-1)); 340 | CHKERRQ (ierr); 341 | ierr = TSSetSolution (thets, x); CHKERRQ (ierr); 342 | 343 | /* Customize timestepping solver, default to Crank-Nicholson */ 344 | ierr = TSSetDuration (thets, time_steps_max, time_total_max); CHKERRQ (ierr); 345 | ierr = TSSetType (thets, TS_CRANK_NICHOLSON); CHKERRQ (ierr); 346 | ierr = TSSetFromOptions (thets); CHKERRQ (ierr); 347 | 348 | /* Run the solver */ 349 | DPRINTF ("Running the solver\n",0); 350 | ierr = TSStep (thets, &steps, &ftime); CHKERRQ (ierr); 351 | 352 | /* Final clean-up */ 353 | DPRINTF ("Cleaning up\n",0); 354 | if (!user.no_contours) 355 | { 356 | if (user.threedee) 357 | { 358 | ierr = GeomviewEnd (PETSC_COMM_WORLD); CHKERRQ (ierr); 359 | } 360 | else 361 | { 362 | ierr = PetscViewerDestroy (user.theviewer); CHKERRQ (ierr); 363 | } 364 | } 365 | ierr = VecDestroy (user.localX); CHKERRQ (ierr); 366 | ierr = VecDestroy (x); CHKERRQ (ierr); 367 | ierr = VecDestroy (user.localF); CHKERRQ (ierr); 368 | ierr = TSDestroy (thets); CHKERRQ (ierr); 369 | ierr = PetscFree (user.label); CHKERRQ (ierr); 370 | ierr = DADestroy (user.da); CHKERRQ (ierr); 371 | 372 | PetscFinalize (); 373 | printf ("Game over man!\n"); 374 | return 0; 375 | } 376 | 377 | 378 | #undef __FUNCT__ 379 | #define __FUNCT__ "FormInitialCondition" 380 | 381 | /*++++++++++++++++++++++++++++++++++++++ 382 | Like it says, put together the initial condition. 383 | 384 | int FormInitialCondition Returns zero or error. 385 | 386 | AppCtx *user The user context structure. 387 | 388 | Vec X Vector in which to place the initial condition. 389 | ++++++++++++++++++++++++++++++++++++++*/ 390 | 391 | int FormInitialCondition (AppCtx *user, Vec X) 392 | { 393 | int i,j,k, row, mc, chvar, mx,my,mz, ierr, xs,ys,zs, xm,ym,zm; 394 | int gxm,gym,gzm, gxs,gys,gzs; 395 | PetscScalar mparam; 396 | PetscScalar *x; 397 | Vec localX = user->localX; 398 | PetscRandom rand; 399 | 400 | mc = user->mc; 401 | chvar = user->chvar; 402 | mx = user->mx; my = user->my; mz = user->mz; 403 | mparam = user->mparam; 404 | 405 | /* Get a pointer to vector data. 406 | - For default PETSc vectors, VecGetArray() returns a pointer to 407 | the data array. Otherwise, the routine is implementation dependent. 408 | - You MUST call VecRestoreArray() when you no longer need access to 409 | the array. */ 410 | if (user->random) 411 | { 412 | DPRINTF ("Setting initial condition as random numbers\n",0); 413 | ierr = PetscRandomCreate (user->comm, RANDOM_DEFAULT, &rand); 414 | CHKERRQ (ierr); 415 | ierr = PetscRandomSetInterval (rand, 0.35, 0.65); CHKERRQ (ierr); 416 | ierr = VecSetRandom (rand, X); CHKERRQ (ierr); 417 | ierr = PetscRandomDestroy (rand); CHKERRQ (ierr); 418 | } 419 | else if (user->load_data > -1) 420 | { 421 | int usermetacount; 422 | char basename [1000], **usermetanames, **usermetadata; 423 | sprintf (basename, "chtsout.time%.3d", user->load_data); 424 | DPRINTF ("Loading data for timestep %d, basename %s\n", user->load_data, 425 | basename); 426 | IlluMultiRead (user->da, X, basename, &usermetacount, &usermetanames, 427 | &usermetadata); 428 | /* TODO: free these */ 429 | for (i=0; i<usermetacount; i++) 430 | PetscPrintf (PETSC_COMM_WORLD, "Parameter \"%s\" = \"%s\"\n", 431 | usermetanames [i], usermetadata [i]); 432 | } 433 | else 434 | { 435 | DPRINTF ("Getting local array\n",0); 436 | ierr = VecGetArray(localX,&x); CHKERRQ (ierr); 437 | 438 | /* Get local grid boundaries (for 2-dimensional DA): 439 | xs, ys - starting grid indices (no ghost points) 440 | xm, ym - widths of local grid (no ghost points) 441 | gxs, gys - starting grid indices (including ghost points) 442 | gxm, gym - widths of local grid (including ghost points) */ 443 | DPRINTF ("Getting corners and ghost corners\n",0); 444 | ierr = DAGetCorners (user->da, &xs,&ys,&zs, &xm,&ym,&zm); CHKERRQ (ierr); 445 | ierr = DAGetGhostCorners (user->da, &gxs,&gys,&gzs, &gxm,&gym,&gzm); 446 | CHKERRQ (ierr); 447 | 448 | /* Set up 2-D so it works */ 449 | if (!user->threedee) { zs = gzs = 0; zm = 1; mz = 10; } 450 | 451 | /* Compute initial condition over the locally owned part of the grid 452 | Initial condition is motionless fluid and equilibrium temperature */ 453 | DPRINTF ("Looping to set initial condition\n",0); 454 | for (k=zs; k<zs+zm; k++) 455 | for (j=ys; j<ys+ym; j++) 456 | for (i=xs; i<xs+xm; i++) 457 | { 458 | row = i - gxs + (j - gys)*gxm + (k-gzs)*gxm*gym; 459 | /* x[C(row)] = (i>=mx*0.43921) ? 1.0 : 0.0; */ 460 | x[C(row)] = ((i<mx/3 || i>2*mx/3) && (j<my/3 || j>2*my/3) && 461 | (k<mz/3 || k>2*mz/3)) ? 1.0 : 0.0; 462 | /* x[C(row)] = (i<mx/2 && j<my/2 && k<mz/2) ? 1.0 : 0.0; */ 463 | /* x[V(row)] = 0.0; 464 | x[Omega(row)] = 0.0; 465 | x[Temp(row)] = (mparam>0)*(PetscScalar)(i)/(PetscScalar)(mx-1); */ 466 | } 467 | 468 | /* Restore vector */ 469 | DPRINTF ("Restoring array to local vector\n",0); 470 | ierr = VecRestoreArray (localX,&x); CHKERRQ (ierr); 471 | 472 | /* Insert values into global vector */ 473 | DPRINTF ("Inserting local vector values into global vector\n",0); 474 | ierr = DALocalToGlobal (user->da,localX,INSERT_VALUES,X); CHKERRQ (ierr); 475 | } 476 | 477 | return 0; 478 | } 479 | 480 | 481 | #undef __FUNCT__ 482 | #define __FUNCT__ "InitializeProblem" 483 | 484 | /*++++++++++++++++++++++++++++++++++++++ 485 | This takes the gory details of initialization out of the way (importing 486 | parameters into the user context, etc.). 487 | 488 | int InitializeProblem Returns zero or error. 489 | 490 | AppCtx *user The user context to fill. 491 | 492 | Vec *xvec Vector into which to put the initial condition. 493 | ++++++++++++++++++++++++++++++++++++++*/ 494 | 495 | int InitializeProblem (AppCtx *user, Vec *xvec) 496 | { 497 | int Nx,Ny,Nz; /* number of processors in x-, y- and z- directions */ 498 | int xs,xm,ys,ym,zs,zm, Nlocal,ierr; 499 | Vec xv; 500 | PetscTruth twodee; 501 | 502 | /* Initialize problem parameters */ 503 | DPRINTF ("Initializing user->parameters\n",0); 504 | user->mx = 20; 505 | user->my = 20; 506 | user->mz = 20; 507 | ierr = PetscOptionsGetInt(PETSC_NULL, "-mx", &user->mx, PETSC_NULL); 508 | CHKERRQ (ierr); 509 | ierr = PetscOptionsGetInt(PETSC_NULL, "-my", &user->my, PETSC_NULL); 510 | CHKERRQ (ierr); 511 | ierr = PetscOptionsGetInt(PETSC_NULL, "-mz", &user->mz, PETSC_NULL); 512 | CHKERRQ (ierr); 513 | /* No. of components in the unknown vector and auxiliary vector */ 514 | user->mc = 1; 515 | /* Problem parameters (kappa, epsilon, gamma, and mparam) */ 516 | user->kappa = 1.0; 517 | ierr = PetscOptionsGetScalar(PETSC_NULL, "-kappa", &user->kappa, PETSC_NULL); 518 | CHKERRQ (ierr); 519 | user->epsilon = 1.0/(user->mx-1); 520 | ierr = PetscOptionsGetScalar (PETSC_NULL, "-epsilon", &user->epsilon, 521 | PETSC_NULL); CHKERRQ (ierr); 522 | user->gamma = 1.0; 523 | ierr = PetscOptionsGetScalar(PETSC_NULL, "-gamma", &user->gamma, PETSC_NULL); 524 | CHKERRQ (ierr); 525 | user->mparam = 0.0; 526 | ierr = PetscOptionsGetScalar (PETSC_NULL, "-mparam", &user->mparam, 527 | PETSC_NULL); CHKERRQ (ierr); 528 | ierr = PetscOptionsHasName (PETSC_NULL, "-twodee", &twodee); 529 | user->threedee = !twodee; 530 | CHKERRQ (ierr); 531 | ierr = PetscOptionsHasName (PETSC_NULL, "-printv", &user->print_vecs); 532 | CHKERRQ (ierr); 533 | ierr = PetscOptionsHasName (PETSC_NULL, "-printg", &user->print_grid); 534 | CHKERRQ (ierr); 535 | ierr = PetscOptionsHasName (PETSC_NULL, "-no_contours", &user->no_contours); 536 | CHKERRQ (ierr); 537 | ierr = PetscOptionsHasName (PETSC_NULL, "-random", &user->random); 538 | CHKERRQ (ierr); 539 | ierr = PetscOptionsHasName (PETSC_NULL, "-save_data", &user->save_data); 540 | CHKERRQ (ierr); 541 | user->load_data = -1; 542 | ierr = PetscOptionsGetInt (PETSC_NULL, "-load_data", &user->load_data, 543 | PETSC_NULL); CHKERRQ (ierr); 544 | 545 | /* Create distributed array (DA) to manage parallel grid and vectors 546 | for principal unknowns (x) and governing residuals (f) 547 | Note the stencil width of 2 for this 4th-order equation. */ 548 | DPRINTF ("Creating distributed array\n",0); 549 | Nx = PETSC_DECIDE; 550 | Ny = PETSC_DECIDE; 551 | Nz = PETSC_DECIDE; 552 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Nx", &Nx, PETSC_NULL);CHKERRQ(ierr); 553 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Ny", &Ny, PETSC_NULL);CHKERRQ(ierr); 554 | ierr = PetscOptionsGetInt (PETSC_NULL, "-Nz", &Nz, PETSC_NULL);CHKERRQ(ierr); 555 | if (user->threedee) 556 | { 557 | DPRINTF ("Three Dee!\n",0); 558 | user->period = DA_XYZPERIODIC; 559 | ierr = DACreate3d (PETSC_COMM_WORLD, user->period, DA_STENCIL_BOX, 560 | user->mx, user->my, user->mz, Nx, Ny, Nz, user->mc, 2, 561 | PETSC_NULL, PETSC_NULL, PETSC_NULL, &user->da); 562 | CHKERRQ (ierr); 563 | } 564 | else 565 | { 566 | user->period = DA_XYPERIODIC; 567 | ierr = DACreate2d (PETSC_COMM_WORLD, user->period, DA_STENCIL_BOX, 568 | user->mx, user->my, Nx, Ny, user->mc, 2, 569 | PETSC_NULL, PETSC_NULL, &user->da); CHKERRQ (ierr); 570 | user->mz = Nz = 1; 571 | } 572 | 573 | /* Extract global vector from DA */ 574 | DPRINTF ("Extracting global vector from DA...\n",0); 575 | ierr = DACreateGlobalVector(user->da,&xv); CHKERRQ (ierr); 576 | 577 | /* Label PDE components */ 578 | DPRINTF ("Labeling PDE components\n",0); 579 | ierr = PetscMalloc (user->mc * sizeof(char*), &(user->label));CHKERRQ (ierr); 580 | user->label[0] = "Concentration (C)"; 581 | /* user->label[1] = "Velocity (V)"; 582 | user->label[2] = "Omega"; 583 | user->label[3] = "Temperature"; */ 584 | ierr = DASetFieldName (user->da, 0, user->label[0]); CHKERRQ(ierr); 585 | 586 | user->x_old = 0; 587 | 588 | /* Get local vector */ 589 | DPRINTF ("Getting local vector\n",0); 590 | ierr = DACreateLocalVector (user->da, &user->localX); CHKERRQ (ierr); 591 | 592 | /* Print grid info */ 593 | if (user->print_grid) 594 | { 595 | DPRINTF ("Printing grid information\n",0); 596 | ierr = DAView(user->da,PETSC_VIEWER_STDOUT_SELF); CHKERRQ (ierr); 597 | ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr); 598 | if (!user->threedee) { 599 | zs = 0; zm = 1; } 600 | ierr = PetscPrintf(PETSC_COMM_WORLD, 601 | "global grid: %d X %d X %d with %d components per" 602 | " node ==> global vector dimension %d\n", 603 | user->mx, user->my, user->mz, user->mc, 604 | user->mc*user->mx*user->my*user->mz); 605 | CHKERRQ (ierr); 606 | fflush(stdout); 607 | ierr = VecGetLocalSize (xv, &Nlocal); CHKERRQ (ierr); 608 | ierr = PetscSynchronizedPrintf 609 | (PETSC_COMM_WORLD,"[%d] local grid %d X %d X %d with %d components" 610 | " per node ==> local vector dimension %d\n", 611 | user->rank, xm, ym, zm, user->mc, Nlocal); CHKERRQ (ierr); 612 | ierr = PetscSynchronizedFlush (PETSC_COMM_WORLD); CHKERRQ (ierr); 613 | } 614 | 615 | /* Compute initial condition */ 616 | DPRINTF ("Forming inital condition\n",0); 617 | ierr = FormInitialCondition (user, xv); CHKERRQ (ierr); 618 | 619 | *xvec = xv; 620 | return 0; 621 | }