<div dir="ltr"><br clear="all"><div>Hi all,</div><div><br></div><div><br></div><div>I took the following petsc example (heat equation dynamics using DMDA) :</div><div><br></div><div><a href="http://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex4.c.html">http://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex4.c.html</a></div><div><br></div><div>that compiles and runs in my system (Kubuntu 64-bit kernel 3.13.0-36-generic with petsc 3.5.2) both serially and in parallel as it should.</div><div><br></div><div>I modified it to solve the quantum Schrodinger equation on a 1d lattice of unit size with a periodic time-dependent RHS.</div><div><br></div><div>I am attaching the diff to this email and pasting it below. A prettier version of the diff is available online at : <a href="http://www.diffnow.com/?report=xhuvj">http://www.diffnow.com/?report=xhuvj</a></div><div><br></div><div>The changed code is running okay when serial, but only shows part of the solution (presumably the part that the root process has) when parallel.</div><div><br></div><div>I tried not to make any major structural changes other than modifying the RHS setting function and hard coding the solver routine (plus, I removed the 'exact solution' routine since my problem doesn't have one), but probably goofed somewhere with the communicators. </div><div><br></div><div>Can someone help me figure it out?</div><div><br></div><div><br></div><div>Regards,</div><div>Analabha</div><div><br></div><div><br></div><div>DIFF pasted (file also attached):</div><div><br></div><div><br></div><div><div>0a1,3</div><div>> //PROBLEMS:</div><div>> //FIX PARALLELIZATION PROBLEM</div><div>> </div><div>11,14c14,24</div><div>< Concepts: TS^time-dependent linear problems</div><div>< Concepts: TS^heat equation</div><div>< Concepts: TS^diffusion equation</div><div>< Processors: n</div><div>---</div><div>> Default problem parameters. These are overwritten by argument vectors</div><div>> */</div><div>> #define DEFAULT_SIZE 100</div><div>> #define DEFAULT_ALPHA 0.3</div><div>> #define DEFAULT_AMPL 3.0</div><div>> #define DEFAULT_W 3.0</div><div>> #define TIME_TOTAL_MAX 100.0<span class="" style="white-space:pre"> </span>/* default max total time */</div><div>> #define TIME_STEPS_MAX 1E5<span class="" style="white-space:pre"> </span>/* default max timesteps */</div><div>> /*</div><div>> Default problem parameters. These are always hard coded and never </div><div>> overwritten</div><div>15a26,29</div><div>> #define INIT_STEP 1E-4</div><div>> #define INIT_TIME 0.0</div><div>> #define ABSERROR 1E-9</div><div>> #define RELERROR 1E-9</div><div>19,31c33,39</div><div>< This program solves the one-dimensional heat equation (also called the</div><div>< diffusion equation),</div><div>< u_t = u_xx,</div><div>< on the domain 0 <= x <= 1, with the boundary conditions</div><div>< u(t,0) = 0, u(t,1) = 0,</div><div>< and the initial condition</div><div>< u(0,x) = sin(6*pi*x) + 3*sin(2*pi*x).</div><div>< This is a linear, second-order, parabolic equation.</div><div>< </div><div>< We discretize the right-hand side using finite differences with</div><div>< uniform grid spacing h:</div><div>< u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2)</div><div>< We then demonstrate time evolution using the various TS methods by</div><div>---</div><div>> This program solves the one-dimensional Schroedinger equation in a </div><div>> tightly bound lattice with a periodic time dependence</div><div>> \partial_t u_m(t) = \sum_m J_{lm}(t) u_m(t),</div><div>> on the domain 0 <= m < L, with the boundary conditions</div><div>> u_m(0) = 0 for all m except m=[L/2], where u_m(0) = 1</div><div>> This is a set of linear, first-order ODEs</div><div>> The time evolution using the various TS methods can be run by</div><div>34,46d41</div><div>< </div><div>< We compare the approximate solution with the exact solution, given by</div><div>< u_exact(x,t) = exp(-36*pi*pi*t) * sin(6*pi*x) +</div><div>< 3*exp(-4*pi*pi*t) * sin(2*pi*x)</div><div>< </div><div>< Notes:</div><div>< This code demonstrates the TS solver interface to two variants of</div><div>< linear problems, u_t = f(u,t), namely</div><div>< - time-dependent f: f(u,t) is a function of t</div><div>< - time-independent f: f(u,t) is simply f(u)</div><div>< </div><div>< The uniprocessor version of this code is ts/examples/tutorials/ex3.c</div><div>< </div><div>59a55</div><div>> </div><div>64d59</div><div>< </div><div>75d69</div><div>< Vec solution;<span class="" style="white-space:pre"> </span>/* global exact solution vector */</div><div>77d70</div><div>< PetscReal h;<span class="" style="white-space:pre"> </span>/* mesh width h = 1/(m-1) */</div><div>79,80c72,77</div><div>< PetscViewer viewer1, viewer2;<span class="" style="white-space:pre"> </span>/* viewers for the solution and error */</div><div>< PetscReal norm_2, norm_max;<span class="" style="white-space:pre"> </span>/* error norms */</div><div>---</div><div>> PetscViewer viewer1;<span class="" style="white-space:pre"> </span>/* viewer for the solution */</div><div>> PetscReal norm_2;<span class="" style="white-space:pre"> </span>/* wavefunction norm */</div><div>> Vec randoms;<span class="" style="white-space:pre"> </span>/* random numbers between -1 and 1 */</div><div>> PetscReal alpha;<span class="" style="white-space:pre"> </span>/* disorder strength */</div><div>> PetscReal h0;<span class="" style="white-space:pre"> </span>/* drive amplitude */</div><div>> PetscReal w;<span class="" style="white-space:pre"> </span>/* drive frequency */</div><div>85a83</div><div>> extern PetscInt kdel (PetscInt, PetscInt);</div><div>87,88c85,86</div><div>< extern PetscErrorCode RHSMatrixHeat (TS, PetscReal, Vec, Mat, Mat, void *);</div><div>< extern PetscErrorCode RHSFunctionHeat (TS, PetscReal, Vec, Vec, void *);</div><div>---</div><div>> extern PetscErrorCode RHSMatrixSchrodinger (TS, PetscReal, Vec, Mat, Mat,</div><div>> <span class="" style="white-space:pre"> </span> void *);</div><div>90d87</div><div>< extern PetscErrorCode ExactSolution (PetscReal, Vec, AppCtx *);</div><div>101,102c98,99</div><div>< PetscReal time_total_max = 1.0;<span class="" style="white-space:pre"> </span>/* default max total time */</div><div>< PetscInt time_steps_max = 100;<span class="" style="white-space:pre"> </span>/* default max timesteps */</div><div>---</div><div>> PetscReal time_total_max = TIME_TOTAL_MAX;<span class="" style="white-space:pre"> </span>/* default max total time */</div><div>> PetscInt time_steps_max = TIME_STEPS_MAX;<span class="" style="white-space:pre"> </span>/* default max timesteps */</div><div>108,109c105,107</div><div>< PetscBool flg;</div><div>< TSProblemType tsproblem = TS_LINEAR;</div><div>---</div><div>> PetscReal alpha;<span class="" style="white-space:pre"> </span>/* Disorder strength */</div><div>> PetscReal h0;<span class="" style="white-space:pre"> </span>/* Drive amplitude */</div><div>> PetscReal w;<span class="" style="white-space:pre"> </span>/* drive frequency */</div><div>119c117</div><div>< m = 60;</div><div>---</div><div>> m = DEFAULT_SIZE;</div><div>121a120,132</div><div>> </div><div>> h0 = DEFAULT_AMPL;</div><div>> ierr = PetscOptionsGetReal (NULL, "-a", &h0, NULL);</div><div>> CHKERRQ (ierr);</div><div>> </div><div>> alpha = DEFAULT_ALPHA;</div><div>> ierr = PetscOptionsGetReal (NULL, "-d", &alpha, NULL);</div><div>> CHKERRQ (ierr);</div><div>> </div><div>> w = DEFAULT_W;</div><div>> ierr = PetscOptionsGetReal (NULL, "-w", &w, NULL);</div><div>> CHKERRQ (ierr);</div><div>> </div><div>123a135</div><div>> </div><div>125,127c137,141</div><div>< appctx.h = 1.0 / (m - 1.0);</div><div>< appctx.norm_2 = 0.0;</div><div>< appctx.norm_max = 0.0;</div><div>---</div><div>> appctx.alpha = alpha;</div><div>> appctx.h0 = h0;</div><div>> appctx.w = w;</div><div>> </div><div>> appctx.norm_2 = 1.0;</div><div>163d176</div><div>< create global work vector for storing exact solution.</div><div>167,168d179</div><div>< ierr = VecDuplicate (u, &appctx.solution);</div><div>< CHKERRQ (ierr);</div><div>171c182</div><div>< Set up displays to show graphs of the solution and error</div><div>---</div><div>> Set up displays to show graphs of the solution </div><div>173,189c184,195</div><div>< </div><div>< ierr =</div><div>< PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "", 80, 380, 400, 160,</div><div>< <span class="" style="white-space:pre"> </span> &appctx.viewer1);</div><div>< CHKERRQ (ierr);</div><div>< ierr = PetscViewerDrawGetDraw (appctx.viewer1, 0, &draw);</div><div>< CHKERRQ (ierr);</div><div>< ierr = PetscDrawSetDoubleBuffer (draw);</div><div>< CHKERRQ (ierr);</div><div>< ierr =</div><div>< PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "", 80, 0, 400, 160,</div><div>< <span class="" style="white-space:pre"> </span> &appctx.viewer2);</div><div>< CHKERRQ (ierr);</div><div>< ierr = PetscViewerDrawGetDraw (appctx.viewer2, 0, &draw);</div><div>< CHKERRQ (ierr);</div><div>< ierr = PetscDrawSetDoubleBuffer (draw);</div><div>< CHKERRQ (ierr);</div><div>---</div><div>> if (appctx.debug)</div><div>> {</div><div>> ierr =</div><div>> <span class="" style="white-space:pre"> </span>PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "Probability Distribution",</div><div>> <span class="" style="white-space:pre"> </span> PETSC_DECIDE, PETSC_DECIDE, PETSC_DRAW_HALF_SIZE,</div><div>> <span class="" style="white-space:pre"> </span> PETSC_DRAW_HALF_SIZE, &appctx.viewer1);</div><div>> CHKERRQ (ierr);</div><div>> ierr = PetscViewerDrawGetDraw (appctx.viewer1, 0, &draw);</div><div>> CHKERRQ (ierr);</div><div>> ierr = PetscDrawSetDoubleBuffer (draw);</div><div>> CHKERRQ (ierr);</div><div>> }</div><div>197,202c203</div><div>< </div><div>< flg = PETSC_FALSE;</div><div>< ierr = PetscOptionsGetBool (NULL, "-nonlinear", &flg, NULL);</div><div>< CHKERRQ (ierr);</div><div>< ierr = TSSetProblemType (ts, flg ? TS_NONLINEAR : TS_LINEAR);</div><div>< CHKERRQ (ierr);</div><div>---</div><div>> TSSetProblemType (ts, TS_LINEAR);<span class="" style="white-space:pre"> </span>/* In this case, the dynamics is always linear */</div><div>205c206</div><div>< Set optional user-defined monitoring routine</div><div>---</div><div>> Set user-defined monitoring routine</div><div>211d211</div><div>< </div><div>224,225c224,257</div><div>< flg = PETSC_FALSE;</div><div>< ierr = PetscOptionsGetBool (NULL, "-time_dependent_rhs", &flg, NULL);</div><div>---</div><div>> /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -</div><div>> Create random number context and generate random numbers </div><div>> into application context</div><div>> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */</div><div>> PetscRandom rctx;<span class="" style="white-space:pre"> </span>/* random number generator context */</div><div>> PetscScalar randomscale = 2.0 * PETSC_i * alpha;</div><div>> </div><div>> ierr = PetscRandomCreate (PETSC_COMM_SELF, &rctx);</div><div>> CHKERRQ (ierr);</div><div>> ierr = PetscRandomSetFromOptions (rctx);</div><div>> CHKERRQ (ierr);</div><div>> ierr = PetscRandomSetInterval (rctx, -1.0, 1.0);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecCreate (PETSC_COMM_SELF, &appctx.randoms);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecSetSizes (appctx.randoms, PETSC_DECIDE, m);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecSetFromOptions (appctx.randoms);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecSetRandom (appctx.randoms, rctx);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecScale (appctx.randoms, randomscale);</div><div>> CHKERRQ (ierr);</div><div>> ierr = PetscRandomDestroy (&rctx);</div><div>> CHKERRQ (ierr);</div><div>> </div><div>> /*</div><div>> For linear problems with a time-dependent f(u,t) in the equation</div><div>> u_t = f(u,t), the user provides the discretized right-hand-side</div><div>> as a time-dependent matrix.</div><div>> */</div><div>> ierr = TSSetRHSFunction (ts, NULL, TSComputeRHSFunctionLinear, &appctx);</div><div>> CHKERRQ (ierr);</div><div>> ierr = TSSetRHSJacobian (ts, A, A, RHSMatrixSchrodinger, &appctx);</div><div>227,254d258</div><div>< if (flg)</div><div>< {</div><div>< /*</div><div>< For linear problems with a time-dependent f(u,t) in the equation</div><div>< u_t = f(u,t), the user provides the discretized right-hand-side</div><div>< as a time-dependent matrix.</div><div>< */</div><div>< ierr = TSSetRHSFunction (ts, NULL, TSComputeRHSFunctionLinear, &appctx);</div><div>< CHKERRQ (ierr);</div><div>< ierr = TSSetRHSJacobian (ts, A, A, RHSMatrixHeat, &appctx);</div><div>< CHKERRQ (ierr);</div><div>< }</div><div>< else</div><div>< {</div><div>< /*</div><div>< For linear problems with a time-independent f(u) in the equation</div><div>< u_t = f(u), the user provides the discretized right-hand-side</div><div>< as a matrix only once, and then sets a null matrix evaluation</div><div>< routine.</div><div>< */</div><div>< ierr = RHSMatrixHeat (ts, 0.0, u, A, A, &appctx);</div><div>< CHKERRQ (ierr);</div><div>< ierr = TSSetRHSFunction (ts, NULL, TSComputeRHSFunctionLinear, &appctx);</div><div>< CHKERRQ (ierr);</div><div>< ierr =</div><div>< <span class="" style="white-space:pre"> </span>TSSetRHSJacobian (ts, A, A, TSComputeRHSJacobianConstant, &appctx);</div><div>< CHKERRQ (ierr);</div><div>< }</div><div>256,266d259</div><div>< if (tsproblem == TS_NONLINEAR)</div><div>< {</div><div>< SNES snes;</div><div>< ierr = TSSetRHSFunction (ts, NULL, RHSFunctionHeat, &appctx);</div><div>< CHKERRQ (ierr);</div><div>< ierr = TSGetSNES (ts, &snes);</div><div>< CHKERRQ (ierr);</div><div>< ierr =</div><div>< <span class="" style="white-space:pre"> </span>SNESSetJacobian (snes, NULL, NULL, SNESComputeJacobianDefault, NULL);</div><div>< CHKERRQ (ierr);</div><div>< }</div><div>272,273c265,266</div><div>< dt = appctx.h * appctx.h / 2.0;</div><div>< ierr = TSSetInitialTimeStep (ts, 0.0, dt);</div><div>---</div><div>> dt = INIT_STEP;</div><div>> ierr = TSSetInitialTimeStep (ts, INIT_TIME, dt);</div><div>281a275</div><div>> - Set default tolerances</div><div>287a282,285</div><div>> ierr = TSSetType (ts, TSRK);</div><div>> CHKERRQ (ierr);</div><div>> ierr = TSRKSetType (ts, TSRK5DP);<span class="" style="white-space:pre"> </span>//Default solver is Runge Kutta (5) Prince-Dormand (4)</div><div>> CHKERRQ (ierr);</div><div>289a288,291</div><div>> </div><div>> ierr = TSSetTolerances (ts, ABSERROR, NULL, RELERROR, NULL);</div><div>> CHKERRQ (ierr);</div><div>> </div><div>321,324c323,324</div><div>< PetscPrintf (PETSC_COMM_WORLD,</div><div>< <span class="" style="white-space:pre"> </span> "Avg. error (2 norm) = %g Avg. error (max norm) = %g\n",</div><div>< <span class="" style="white-space:pre"> </span> (double) (appctx.norm_2 / steps),</div><div>< <span class="" style="white-space:pre"> </span> (double) (appctx.norm_max / steps));</div><div>---</div><div>> PetscPrintf (PETSC_COMM_WORLD, "Final. (2 norm) = %g \n",</div><div>> <span class="" style="white-space:pre"> </span> (double) (appctx.norm_2));</div><div>338,341c338,343</div><div>< ierr = PetscViewerDestroy (&appctx.viewer1);</div><div>< CHKERRQ (ierr);</div><div>< ierr = PetscViewerDestroy (&appctx.viewer2);</div><div>< CHKERRQ (ierr);</div><div>---</div><div>> </div><div>> if (appctx.debug)</div><div>> {</div><div>> ierr = PetscViewerDestroy (&appctx.viewer1);</div><div>> CHKERRQ (ierr);</div><div>> }</div><div>344,345d345</div><div>< ierr = VecDestroy (&appctx.solution);</div><div>< CHKERRQ (ierr);</div><div>347a348,349</div><div>> ierr = VecDestroy (&appctx.randoms);</div><div>> CHKERRQ (ierr);</div><div>377c379</div><div>< PetscScalar *u_localptr, h = appctx->h;</div><div>---</div><div>> PetscScalar *u_localptr;</div><div>379a382</div><div>> PetscScalar zero = 0.0, one = 1.0;</div><div>394,395d396</div><div>< - Note that the Fortran interface to VecGetArray() differs from the</div><div>< C version. See the users manual for details.</div><div>402,403c403,404</div><div>< directly into the array locations. Alternatively, we could use</div><div>< VecSetValues() or VecSetValuesLocal().</div><div>---</div><div>> directly into the array locations.</div><div>> This will set the midpoint to unity</div><div>406,408c407,410</div><div>< u_localptr[i - mybase] =</div><div>< PetscSinScalar (PETSC_PI * i * 6. * h) +</div><div>< 3. * PetscSinScalar (PETSC_PI * i * 2. * h);</div><div>---</div><div>> if (i - mybase == floor ((appctx->m) / 2))</div><div>> u_localptr[i - mybase] = one;</div><div>> else</div><div>> u_localptr[i - mybase] = zero;</div><div>421c423</div><div>< ierr = PetscPrintf (appctx->comm, "initial guess vector\n");</div><div>---</div><div>> ierr = PetscPrintf (appctx->comm, "initial vector\n");</div><div>432,486d433</div><div>< #define __FUNCT__ "ExactSolution"</div><div>< /*</div><div>< ExactSolution - Computes the exact solution at a given time.</div><div>< </div><div>< Input Parameters:</div><div>< t - current time</div><div>< solution - vector in which exact solution will be computed</div><div>< appctx - user-defined application context</div><div>< </div><div>< Output Parameter:</div><div>< solution - vector with the newly computed exact solution</div><div>< */</div><div>< PetscErrorCode</div><div>< ExactSolution (PetscReal t, Vec solution, AppCtx * appctx)</div><div>< {</div><div>< PetscScalar *s_localptr, h = appctx->h, ex1, ex2, sc1, sc2;</div><div>< PetscInt i, mybase, myend;</div><div>< PetscErrorCode ierr;</div><div>< </div><div>< /*</div><div>< Determine starting and ending points of each processor's</div><div>< range of grid values</div><div>< */</div><div>< ierr = VecGetOwnershipRange (solution, &mybase, &myend);</div><div>< CHKERRQ (ierr);</div><div>< </div><div>< /*</div><div>< Get a pointer to vector data.</div><div>< */</div><div>< ierr = VecGetArray (solution, &s_localptr);</div><div>< CHKERRQ (ierr);</div><div>< </div><div>< /*</div><div>< Simply write the solution directly into the array locations.</div><div>< Alternatively, we culd use VecSetValues() or VecSetValuesLocal().</div><div>< */</div><div>< ex1 = PetscExpReal (-36. * PETSC_PI * PETSC_PI * t);</div><div>< ex2 = PetscExpReal (-4. * PETSC_PI * PETSC_PI * t);</div><div>< sc1 = PETSC_PI * 6. * h;</div><div>< sc2 = PETSC_PI * 2. * h;</div><div>< for (i = mybase; i < myend; i++)</div><div>< s_localptr[i - mybase] =</div><div>< PetscSinScalar (sc1 * (PetscReal) i) * ex1 +</div><div>< 3. * PetscSinScalar (sc2 * (PetscReal) i) * ex2;</div><div>< </div><div>< /*</div><div>< Restore vector</div><div>< */</div><div>< ierr = VecRestoreArray (solution, &s_localptr);</div><div>< CHKERRQ (ierr);</div><div>< return 0;</div><div>< }</div><div>< </div><div>< /* --------------------------------------------------------------------- */</div><div>< #undef __FUNCT__</div><div>490,491c437</div><div>< each timestep. This example plots the solution and computes the</div><div>< error in two different norms.</div><div>---</div><div>> each timestep. This example plots the solution and returns the norm</div><div>509c455,457</div><div>< PetscReal norm_2, norm_max;</div><div>---</div><div>> PetscReal norm_2;</div><div>> PetscScalar edge_mobility, variance;</div><div>> PetscInt ix = (appctx->m) - 1;</div><div>512c460</div><div>< View a graph of the current iterate</div><div>---</div><div>> Compute and store norm of wavefunction in norm_2 </div><div>514c462</div><div>< ierr = VecView (u, appctx->viewer2);</div><div>---</div><div>> ierr = VecNorm (u, NORM_2, &norm_2);</div><div>515a464</div><div>> appctx->norm_2 = norm_2;</div><div>518c467</div><div>< Compute the exact solution</div><div>---</div><div>> Compute variance</div><div>520c469,476</div><div>< ierr = ExactSolution (time, appctx->solution, appctx);</div><div>---</div><div>> </div><div>> /* Build the lattice sites */</div><div>> Vec uabs, sites, sitesq;</div><div>> ierr = VecDuplicate (u, &sites);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecDuplicate (u, &sitesq);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecDuplicate (u, &uabs);</div><div>523,526c479,483</div><div>< /*</div><div>< Print debugging information if desired</div><div>< */</div><div>< if (appctx->debug)</div><div>---</div><div>> PetscInt c;</div><div>> PetscScalar cdbl;</div><div>> PetscScalar csq;</div><div>> </div><div>> for (c = 1; c < appctx->m; c++)</div><div>528c485,486</div><div>< ierr = PetscPrintf (appctx->comm, "Computed solution vector\n");</div><div>---</div><div>> cdbl = c;</div><div>> ierr = VecSetValues (sites, 1, &c, &cdbl, INSERT_VALUES);</div><div>530,534c488,489</div><div>< ierr = VecView (u, PETSC_VIEWER_STDOUT_WORLD);</div><div>< CHKERRQ (ierr);</div><div>< ierr = PetscPrintf (appctx->comm, "Exact solution vector\n");</div><div>< CHKERRQ (ierr);</div><div>< ierr = VecView (appctx->solution, PETSC_VIEWER_STDOUT_WORLD);</div><div>---</div><div>> csq = c * c;</div><div>> ierr = VecSetValues (sitesq, 1, &c, &csq, INSERT_VALUES);</div><div>538,541c493,495</div><div>< /*</div><div>< Compute the 2-norm and max-norm of the error</div><div>< */</div><div>< ierr = VecAXPY (appctx->solution, -1.0, u);</div><div>---</div><div>> ierr = VecAssemblyBegin (sites);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecAssemblyEnd (sites);</div><div>543c497</div><div>< ierr = VecNorm (appctx->solution, NORM_2, &norm_2);</div><div>---</div><div>> ierr = VecAssemblyBegin (sitesq);</div><div>545,546c499</div><div>< norm_2 = PetscSqrtReal (appctx->h) * norm_2;</div><div>< ierr = VecNorm (appctx->solution, NORM_MAX, &norm_max);</div><div>---</div><div>> ierr = VecAssemblyEnd (sitesq);</div><div>549,556c502,507</div><div>< /*</div><div>< PetscPrintf() causes only the first processor in this</div><div>< communicator to print the timestep information.</div><div>< */</div><div>< ierr =</div><div>< PetscPrintf (appctx->comm,</div><div>< <span class="" style="white-space:pre"> </span> "Timestep %D: time = %g 2-norm error = %g max norm error = %g\n",</div><div>< <span class="" style="white-space:pre"> </span> step, (double) time, (double) norm_2, (double) norm_max);</div><div>---</div><div>> ierr = VecCopy (u, uabs);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecAbs (uabs);</div><div>> CHKERRQ (ierr);</div><div>> </div><div>> ierr = VecPointwiseMult (uabs, uabs, uabs);</div><div>558,559c509,520</div><div>< appctx->norm_2 += norm_2;</div><div>< appctx->norm_max += norm_max;</div><div>---</div><div>> </div><div>> /*Compute \sum_m m^2 |u_m|^2 */</div><div>> PetscScalar msqbar;</div><div>> ierr = VecTDot (uabs, sitesq, &msqbar);</div><div>> CHKERRQ (ierr);</div><div>> /*Compute (\sum_m m |u_m|^2)^2 */</div><div>> PetscScalar mbarsq;</div><div>> ierr = VecTDot (uabs, sites, &mbarsq);</div><div>> CHKERRQ (ierr);</div><div>> mbarsq = mbarsq * mbarsq;</div><div>> variance = msqbar - mbarsq;</div><div>> variance = variance / ((appctx->m) * (appctx->m));</div><div>562c523</div><div>< View a graph of the error</div><div>---</div><div>> Compute Edge Mobility</div><div>564c525</div><div>< ierr = VecView (appctx->solution, appctx->viewer1);</div><div>---</div><div>> ierr = VecGetValues (u, 1, &ix, &edge_mobility);</div><div>565a527,554</div><div>> PetscReal edge_mob_ampl = cabs (edge_mobility);</div><div>> </div><div>> if (time == INIT_TIME && !appctx->debug)</div><div>> {</div><div>> ierr =</div><div>> <span class="" style="white-space:pre"> </span>PetscPrintf (PETSC_COMM_WORLD,</div><div>> <span class="" style="white-space:pre"> </span> " ###########################################\n", time,</div><div>> <span class="" style="white-space:pre"> </span> edge_mob_ampl);</div><div>> CHKERRQ (ierr);</div><div>> ierr =</div><div>> <span class="" style="white-space:pre"> </span>PetscPrintf (PETSC_COMM_WORLD,</div><div>> <span class="" style="white-space:pre"> </span> " # time \t# edge_mob\t# variance #\n", time,</div><div>> <span class="" style="white-space:pre"> </span> edge_mob_ampl);</div><div>> CHKERRQ (ierr);</div><div>> ierr =</div><div>> <span class="" style="white-space:pre"> </span>PetscPrintf (PETSC_COMM_WORLD,</div><div>> <span class="" style="white-space:pre"> </span> " ###########################################\n", time,</div><div>> <span class="" style="white-space:pre"> </span> edge_mob_ampl);</div><div>> CHKERRQ (ierr);</div><div>> </div><div>> }</div><div>> if (!appctx->debug)</div><div>> {</div><div>> ierr =</div><div>> <span class="" style="white-space:pre"> </span>PetscPrintf (PETSC_COMM_WORLD, " # %3.5lf\t# %.5lf\t# %3.5lf #\n",</div><div>> <span class="" style="white-space:pre"> </span> time, edge_mob_ampl, variance);</div><div>> CHKERRQ (ierr);</div><div>> }</div><div>572,574c561,568</div><div>< ierr = PetscPrintf (appctx->comm, "Error vector\n");</div><div>< CHKERRQ (ierr);</div><div>< ierr = VecView (appctx->solution, PETSC_VIEWER_STDOUT_WORLD);</div><div>---</div><div>> /*</div><div>> View a graph of the current iterate</div><div>> */</div><div>> ierr =</div><div>> <span class="" style="white-space:pre"> </span>PetscPrintf (PETSC_COMM_WORLD, "time = %lf\t norm = %lf\n", time,</div><div>> <span class="" style="white-space:pre"> </span> norm_2);</div><div>> </div><div>> ierr = VecView (uabs, appctx->viewer1);</div><div>576a571,576</div><div>> ierr = VecDestroy (&sites);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecDestroy (&sitesq);</div><div>> CHKERRQ (ierr);</div><div>> ierr = VecDestroy (&uabs);</div><div>> CHKERRQ (ierr);</div><div>583c583</div><div>< #define __FUNCT__ "RHSMatrixHeat"</div><div>---</div><div>> #define __FUNCT__ "RHSMatrixSchrodinger"</div><div>585c585</div><div>< RHSMatrixHeat - User-provided routine to compute the right-hand-side</div><div>---</div><div>> RHSMatrixSchrodinger - User-provided routine to compute the right-hand-side</div><div>599,600c599</div><div>< Notes:</div><div>< RHSMatrixHeat computes entries for the locally owned part of the system.</div><div>---</div><div>> RHSMatrixSchrodinger computes entries for the locally owned part of the system.</div><div>613c612</div><div>< RHSMatrixHeat (TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)</div><div>---</div><div>> RHSMatrixSchrodinger (TS ts, PetscReal t, Vec X, Mat AA, Mat BB, void *ctx)</div><div>618,619c617,629</div><div>< PetscInt i, mstart, mend, idx[3];</div><div>< PetscScalar v[3], stwo = -2. / (appctx->h * appctx->h), sone = -.5 * stwo;</div><div>---</div><div>> PetscInt m = appctx->m;</div><div>> PetscScalar drive;</div><div>> PetscReal h0 = appctx->h0;</div><div>> PetscReal w = appctx->w;</div><div>> PetscScalar *randoms_loc;</div><div>> </div><div>> ierr = VecGetArray (appctx->randoms, &randoms_loc);</div><div>> CHKERRQ (ierr);</div><div>> /*Calculate drive amplitude at this time */</div><div>> </div><div>> drive = w * t;</div><div>> drive = PetscSinReal (drive);</div><div>> drive = 2.0 * PETSC_i * h0 * drive;</div><div>623,625c633,641</div><div>< - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */</div><div>< </div><div>< ierr = MatGetOwnershipRange (A, &mstart, &mend);</div><div>---</div><div>> The Jacobian Matrix in our case is blocked</div><div>> - The values are A_{lm} = i\left(\delta_{lm+1}+\delta_{lm-1}\right)+2i \alpha h_m\delta_{lm} + 2i h(t) \delta_{lm}</div><div>> - Thus the superdiagonal and subdiagonal elements are just i, and diagonal elements are random + drive</div><div>> - The elements are set rowwise.</div><div>> - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */</div><div>> PetscInt rowstart, rowend, rowcount;</div><div>> PetscScalar edge[2], middle[3];</div><div>> PetscInt idxedge[2], idxmiddle[3];</div><div>> ierr = MatGetOwnershipRange (A, &rowstart, &rowend);</div><div>627,655c643</div><div>< </div><div>< /*</div><div>< Set matrix rows corresponding to boundary data</div><div>< */</div><div>< </div><div>< if (mstart == 0)</div><div>< {<span class="" style="white-space:pre"> </span>/* first processor only */</div><div>< v[0] = 1.0;</div><div>< ierr = MatSetValues (A, 1, &mstart, 1, &mstart, v, INSERT_VALUES);</div><div>< CHKERRQ (ierr);</div><div>< mstart++;</div><div>< }</div><div>< </div><div>< if (mend == appctx->m)</div><div>< {<span class="" style="white-space:pre"> </span>/* last processor only */</div><div>< mend--;</div><div>< v[0] = 1.0;</div><div>< ierr = MatSetValues (A, 1, &mend, 1, &mend, v, INSERT_VALUES);</div><div>< CHKERRQ (ierr);</div><div>< }</div><div>< </div><div>< /*</div><div>< Set matrix rows corresponding to interior data. We construct the</div><div>< matrix one row at a time.</div><div>< */</div><div>< v[0] = sone;</div><div>< v[1] = stwo;</div><div>< v[2] = sone;</div><div>< for (i = mstart; i < mend; i++)</div><div>---</div><div>> for (rowcount = rowstart; rowcount < rowend; rowcount++)</div><div>657,661c645,677</div><div>< idx[0] = i - 1;</div><div>< idx[1] = i;</div><div>< idx[2] = i + 1;</div><div>< ierr = MatSetValues (A, 1, &i, 3, idx, v, INSERT_VALUES);</div><div>< CHKERRQ (ierr);</div><div>---</div><div>> if (rowcount == 0)<span class="" style="white-space:pre"> </span>//If it is the first row of the whole matrix</div><div>> <span class="" style="white-space:pre"> </span>{</div><div>> <span class="" style="white-space:pre"> </span> edge[0] = drive + randoms_loc[0];</div><div>> <span class="" style="white-space:pre"> </span> edge[1] = PETSC_i;</div><div>> <span class="" style="white-space:pre"> </span> idxedge[0] = 0;</div><div>> <span class="" style="white-space:pre"> </span> idxedge[1] = 1;</div><div>> <span class="" style="white-space:pre"> </span> ierr =</div><div>> <span class="" style="white-space:pre"> </span> MatSetValues (A, 1, &rowcount, 2, idxedge, edge, INSERT_VALUES);</div><div>> <span class="" style="white-space:pre"> </span> CHKERRQ (ierr);</div><div>> <span class="" style="white-space:pre"> </span>}</div><div>> else if (rowcount == m - 1)<span class="" style="white-space:pre"> </span>//If it is the last row of the whole matrix</div><div>> <span class="" style="white-space:pre"> </span>{</div><div>> <span class="" style="white-space:pre"> </span> edge[0] = PETSC_i;</div><div>> <span class="" style="white-space:pre"> </span> edge[1] = drive + randoms_loc[m - 1];</div><div>> <span class="" style="white-space:pre"> </span> idxedge[0] = m - 2;</div><div>> <span class="" style="white-space:pre"> </span> idxedge[1] = m - 1;</div><div>> <span class="" style="white-space:pre"> </span> ierr =</div><div>> <span class="" style="white-space:pre"> </span> MatSetValues (A, 1, &rowcount, 2, idxedge, edge, INSERT_VALUES);</div><div>> <span class="" style="white-space:pre"> </span> CHKERRQ (ierr);</div><div>> <span class="" style="white-space:pre"> </span>}</div><div>> else</div><div>> <span class="" style="white-space:pre"> </span>{</div><div>> <span class="" style="white-space:pre"> </span> middle[0] = PETSC_i;</div><div>> <span class="" style="white-space:pre"> </span> middle[1] = drive + randoms_loc[rowcount];</div><div>> <span class="" style="white-space:pre"> </span> middle[2] = PETSC_i;</div><div>> <span class="" style="white-space:pre"> </span> idxmiddle[0] = rowcount - 1;</div><div>> <span class="" style="white-space:pre"> </span> idxmiddle[1] = rowcount;</div><div>> <span class="" style="white-space:pre"> </span> idxmiddle[2] = rowcount + 1;</div><div>> <span class="" style="white-space:pre"> </span> ierr =</div><div>> <span class="" style="white-space:pre"> </span> MatSetValues (A, 1, &rowcount, 3, idxmiddle, middle,</div><div>> <span class="" style="white-space:pre"> </span> INSERT_VALUES);</div><div>> <span class="" style="white-space:pre"> </span> CHKERRQ (ierr);</div><div>> <span class="" style="white-space:pre"> </span>}</div><div>683a700,701</div><div>> ierr = VecRestoreArray (appctx->randoms, &randoms_loc);</div><div>> CHKERRQ (ierr);</div><div>689,691c707,709</div><div>< #define __FUNCT__ "RHSFunctionHeat"</div><div>< PetscErrorCode</div><div>< RHSFunctionHeat (TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)</div><div>---</div><div>> #define __FUNCT__ "kdel"</div><div>> PetscInt</div><div>> kdel (PetscInt i, PetscInt j)</div><div>693,704c711,714</div><div>< PetscErrorCode ierr;</div><div>< Mat A;</div><div>< </div><div>< PetscFunctionBeginUser;</div><div>< ierr = TSGetRHSJacobian (ts, &A, NULL, NULL, &ctx);</div><div>< CHKERRQ (ierr);</div><div>< ierr = RHSMatrixHeat (ts, t, globalin, A, NULL, ctx);</div><div>< CHKERRQ (ierr);</div><div>< /* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */</div><div>< ierr = MatMult (A, globalin, globalout);</div><div>< CHKERRQ (ierr);</div><div>< PetscFunctionReturn (0);</div><div>---</div><div>> PetscInt value = 0;</div><div>> if (i == j)</div><div>> value = 1;</div><div>> return value;</div></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div>-- <br><font face="tahoma, sans-serif">---</font><div><font face="tahoma, sans-serif"><b>Analabha Roy</b></font></div><div><font face="tahoma, sans-serif"><a href="http://www.csir.res.in" target="_blank">C.S.I.R</a> <a href="http://csirhrdg.res.in/poolsra.htm" target="_blank">Senior Research Associate</a></font></div><div><font face="tahoma, sans-serif"><a href="http://www.saha.ac.in" target="_blank">Saha Institute of Nuclear Physics</a></font></div><div><font face="tahoma, sans-serif">Section 1, Block AF</font></div><div><font face="tahoma, sans-serif">Bidhannagar, Calcutta 700064</font></div><div><font face="tahoma, sans-serif">India</font></div><div><b style="font-family:tahoma,sans-serif">Emails</b><span style="font-family:tahoma,sans-serif">: </span><a href="mailto:daneel@physics.utexas.edu" style="font-family:tahoma,sans-serif" target="_blank">daneel@physics.utexas.edu</a><span style="font-family:tahoma,sans-serif">, </span><a href="mailto:hariseldon99@gmail.com" style="font-family:tahoma,sans-serif" target="_blank">hariseldon99@gmail.com</a></div><div><font face="tahoma, sans-serif"><b>Webpage</b>: <a href="http://www.ph.utexas.edu/~daneel/" target="_blank">http://www.ph.utexas.edu/~daneel/</a></font></div>
</div>