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