<div dir="ltr">Hi all,<div><br></div><div> Update: Turns out i set my initial condition wrong due to incorrectly using local indices while setting the value, and running it on multiple processor causes overstriding, setting all vector elements to zero. </div><div><br></div><div>I have fixed it. Basically, line 407 changes from the local index (i-ibase) to the global index (just i).</div><div><br></div><div>  for (i = mybase; i < myend; i++)</div><div>    if (i == floor ((appctx->m) / 2))</div><div>      u_localptr[i - mybase] = one;</div><div>    else</div><div>      u_localptr[i - mybase] = zero; </div><div><br></div><div>A trivially silly thing I'm afraid.<br></div><div><br></div><div>Thanks for pointing out my other unrelated error. I think I will compute the variance using MatMult() instead so as to keep everything collective.</div><div><br></div><div>Basically set a diagonal matrix to X_{ii} = i and Xsq_{ii} = i^2 and compute expectation values using MatMult() and VecDot()</div><div><br></div><div>Regards,</div><div>Analabha</div><div><br></div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Oct 16, 2014 at 1:18 PM, Analabha Roy <span dir="ltr"><<a href="mailto:hariseldon99@gmail.com" target="_blank">hariseldon99@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi all,<div><br></div><div>Thanks for the response.</div><div class="gmail_extra"><br><div class="gmail_quote"><span class="">On Wed, Oct 15, 2014 at 4:15 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>I am not sure what you are doing at line 483. You seem to have all processes indexing a Vec</div><div>at all entries, which cannot be correct.</div><div><br></div></div></div></div></blockquote><div><br></div></span><div>Yes, you're right, of course. I have removed that bit (I'll do it in post-processing after I have dumped the full state data to file).</div><div><br></div><div>In the meantime, the graphical output is still messed up when the code is run with debug flag in parallel. However, parallel graphical output works fine in the original ex4.c. I can;t figure out why...</div><div><br></div><div>I have updated diff (<a href="http://www.diffnow.com/?report=49thu" target="_blank">link to pretty version</a>) attached and posted below.</div><div><br></div><div><br></div><div>Thanks and Regards,</div><div>Analabha</div><div><br></div><div><div>0a1,2</div><div>> //PROBLEMS:</div><div>> //FIX PARALLELIZATION PROBLEM</div><div>3,8c5,11</div><div><   "Solves a simple time-dependent linear PDE (the heat equation).\n\</div><div>< Input parameters include:\n\</div><div><   -m <points>, where <points> = number of grid points\n\</div><div><   -time_dependent_rhs : Treat the problem as having a time-dependent right-hand side\n\</div><div><   -debug              : Activate debugging printouts\n\</div><div><   -nox                : Deactivate x-window graphics\n\n";</div><div>---</div><div>>   "Solves a simple time-dependent Schroedinger equation for a TB model with periodic drive.\n\</div><div>> Apart from PetSc's routine options, input parameters include:\n\</div><div>>   -m <pnts>,   where <pnts>   = lattice size\n\</div><div>>   -a <ampl>,   where <ampl>   = drive amplitude\n\</div><div>>   -d <strg>,   where <strg>   = disorder strength\n\</div><div>>   -w <freq>,   where <freq>   = drive frequency\n\</div><div>>   -debug              : Activate debugging printouts\n\n";</div><div>11,14c14</div><span class=""><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></span><div>15a16,29</div><span class=""><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 style="white-space:pre-wrap">      </span>/* default max total time */</div><div>> #define TIME_STEPS_MAX 1E5<span style="white-space:pre-wrap">      </span>/* default max timesteps */</div><div>> /*</div><div>>    Default problem parameters. These are always hard coded and never </div><div>>    overwritten</div></span><div>> */</div><div><div class="h5"><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 style="white-space:pre-wrap">                        </span>/* global exact solution vector */</div><div>77d70</div><div><   PetscReal h;<span style="white-space:pre-wrap">                       </span>/* mesh width h = 1/(m-1) */</div><div>79,80c72,77</div><div><   PetscViewer viewer1, viewer2;<span style="white-space:pre-wrap">      </span>/* viewers for the solution and error */</div><div><   PetscReal norm_2, norm_max;<span style="white-space:pre-wrap">      </span>/* error norms */</div><div>---</div><div>>   PetscViewer viewer1;<span style="white-space:pre-wrap">          </span>/* viewer for the solution  */</div><div>>   PetscReal norm_2;<span style="white-space:pre-wrap">         </span>/* wavefunction norm */</div><div>>   Vec randoms;<span style="white-space:pre-wrap">                      </span>/* random numbers between -1 and 1 */</div><div>>   PetscReal alpha;<span style="white-space:pre-wrap">            </span>/* disorder strength */</div><div>>   PetscReal h0;<span style="white-space:pre-wrap">                     </span>/* drive amplitude */</div><div>>   PetscReal w;<span style="white-space:pre-wrap">                        </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 style="white-space:pre-wrap">                                 </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 style="white-space:pre-wrap"> </span>/* default max total time */</div><div><   PetscInt time_steps_max = 100;<span style="white-space:pre-wrap">       </span>/* default max timesteps */</div><div>---</div><div>>   PetscReal time_total_max = TIME_TOTAL_MAX;<span style="white-space:pre-wrap">  </span>/* default max total time */</div><div>>   PetscInt time_steps_max = TIME_STEPS_MAX;<span style="white-space:pre-wrap">    </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 style="white-space:pre-wrap">              </span>/* Disorder strength */</div><div>>   PetscReal h0;<span style="white-space:pre-wrap">                     </span>/* Drive amplitude */</div><div>>   PetscReal w;<span style="white-space:pre-wrap">                        </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 style="white-space:pre-wrap">                 </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 style="white-space:pre-wrap">                   </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 style="white-space:pre-wrap">  </span>PetscViewerDrawOpen (PETSC_COMM_WORLD, 0, "Probability Distribution",</div><div>> <span style="white-space:pre-wrap">                     </span>     PETSC_DECIDE, PETSC_DECIDE, PETSC_DRAW_HALF_SIZE,</div><div>> <span style="white-space:pre-wrap">                    </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 style="white-space:pre-wrap">      </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 style="white-space:pre-wrap">         </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 style="white-space:pre-wrap">   </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 style="white-space:pre-wrap">      </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 style="white-space:pre-wrap">   </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 style="white-space:pre-wrap">             </span> "Avg. error (2 norm) = %g Avg. error (max norm) = %g\n",</div><div>< <span style="white-space:pre-wrap">         </span> (double) (appctx.norm_2 / steps),</div><div>< <span style="white-space:pre-wrap">          </span> (double) (appctx.norm_max / steps));</div><div>---</div><div>>     PetscPrintf (PETSC_COMM_WORLD, "Final. (2 norm) = %g \n",</div><div>> <span style="white-space:pre-wrap">                </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></div><div>509,515c455</div><div><   PetscReal norm_2, norm_max;</div><span class=""><div>< </div><div><   /*</div><div><      View a graph of the current iterate</div><div><    */</div></span><div><   ierr = VecView (u, appctx->viewer2);</div><div><   CHKERRQ (ierr);</div><div>---</div><div>>   PetscReal norm_2;</div><div>518c458</div><span class=""><div><      Compute the exact solution</div><div>---</div></span><span class=""><div>>      Compute and store norm of wavefunction in norm_2 </div></span><div>520c460</div><span class=""><div><   ierr = ExactSolution (time, appctx->solution, appctx);</div><div>---</div></span><span class=""><div>>   ierr = VecNorm (u, NORM_2, &norm_2);</div></span><div>521a462</div><div>>   appctx->norm_2 = norm_2;</div><div>523,536d463</div><span class=""><div><   /*</div><div><      Print debugging information if desired</div><div><    */</div><div><   if (appctx->debug)</div></span><span class=""><div><     {</div><div><       ierr = PetscPrintf (appctx->comm, "Computed solution vector\n");</div></span><span class=""><div><       CHKERRQ (ierr);</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></span><div><       CHKERRQ (ierr);</div><div><     }</div><div>538,541c465,467</div><span class=""><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></span><span class=""><div>>   /* Build the lattice sites */</div></span><div>>   Vec uabs;</div><span class=""><div>>   ierr = VecDuplicate (u, &uabs);</div></span><div>543c469,470</div><span class=""><div><   ierr = VecNorm (appctx->solution, NORM_2, &norm_2);</div><div>---</div><div>> </div></span><div>>   ierr = VecCopy (u, uabs);</div><div>545,546c472</div><span class=""><div><   norm_2 = PetscSqrtReal (appctx->h) * norm_2;</div><div><   ierr = VecNorm (appctx->solution, NORM_MAX, &norm_max);</div><div>---</div></span><div>>   ierr = VecAbs (uabs);</div><div>549,556c475</div><span class=""><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 style="white-space:pre-wrap">            </span> "Timestep %D: time = %g 2-norm error = %g max norm error = %g\n",</div><div>< <span style="white-space:pre-wrap">                </span> step, (double) time, (double) norm_2, (double) norm_max);</div><div>---</div></span><div>>   ierr = VecPointwiseMult (uabs, uabs, uabs);</div><div>558,559d476</div><span class=""><div><   appctx->norm_2 += norm_2;</div><div><   appctx->norm_max += norm_max;</div></span><div>561,565c478,483</div><div><   /*</div><div><      View a graph of the error</div><div><    */</div><span class=""><div><   ierr = VecView (appctx->solution, appctx->viewer1);</div></span><div><   CHKERRQ (ierr);</div><div>---</div><div>>   //if (!appctx->debug)</div><div>>   //{</div><div>>   /* Create another viewer in application context and dump full output to it */</div><div>>   //ierr = VecView (u, appctx->viewer2);</div><div>>   //CHKERRQ (ierr);</div><div>>   //}</div><div>572,574c490,497</div><span class=""><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 style="white-space:pre-wrap"> </span>PetscPrintf (PETSC_COMM_WORLD, "time = %lf\t norm = %lf\n", time,</div><div>> <span style="white-space:pre-wrap">         </span>     norm_2);</div><div>> </div><div>>       ierr = VecView (uabs, appctx->viewer1);</div></span><div>576a500,501</div><span class=""><div>>   ierr = VecDestroy (&uabs);</div><div>>   CHKERRQ (ierr);</div></span><div>583c508</div><span class=""><div>< #define __FUNCT__ "RHSMatrixHeat"</div><div>---</div><div>> #define __FUNCT__ "RHSMatrixSchrodinger"</div></span><div>585c510</div><span class=""><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></span><div>599,600c524</div><span class=""><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></span><div>613c537</div><span class=""><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></span><div>618,619c542,554</div><span class=""><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></span><div>623,625c558,566</div><span class=""><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></span><div>627,655c568</div><div><div class="h5"><div>< </div><div><   /*</div><div><      Set matrix rows corresponding to boundary data</div><div><    */</div><div>< </div><div><   if (mstart == 0)</div><div><     {<span style="white-space:pre-wrap">                                </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 style="white-space:pre-wrap">                          </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></div><div>657,661c570,602</div><div><div class="h5"><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 style="white-space:pre-wrap">       </span>//If it is the first row of the whole matrix</div><div>> <span style="white-space:pre-wrap">        </span>{</div><div>> <span style="white-space:pre-wrap">   </span>  edge[0] = drive + randoms_loc[0];</div><div>> <span style="white-space:pre-wrap">        </span>  edge[1] = PETSC_i;</div><div>> <span style="white-space:pre-wrap">       </span>  idxedge[0] = 0;</div><div>> <span style="white-space:pre-wrap">  </span>  idxedge[1] = 1;</div><div>> <span style="white-space:pre-wrap">  </span>  ierr =</div><div>> <span style="white-space:pre-wrap">   </span>    MatSetValues (A, 1, &rowcount, 2, idxedge, edge, INSERT_VALUES);</div><div>> <span style="white-space:pre-wrap">  </span>  CHKERRQ (ierr);</div><div>> <span style="white-space:pre-wrap">  </span>}</div><div>>       else if (rowcount == m - 1)<span style="white-space:pre-wrap">       </span>//If it is the last row of the whole matrix</div><div>> <span style="white-space:pre-wrap"> </span>{</div><div>> <span style="white-space:pre-wrap">   </span>  edge[0] = PETSC_i;</div><div>> <span style="white-space:pre-wrap">       </span>  edge[1] = drive + randoms_loc[m - 1];</div><div>> <span style="white-space:pre-wrap">    </span>  idxedge[0] = m - 2;</div><div>> <span style="white-space:pre-wrap">      </span>  idxedge[1] = m - 1;</div><div>> <span style="white-space:pre-wrap">      </span>  ierr =</div><div>> <span style="white-space:pre-wrap">   </span>    MatSetValues (A, 1, &rowcount, 2, idxedge, edge, INSERT_VALUES);</div><div>> <span style="white-space:pre-wrap">  </span>  CHKERRQ (ierr);</div><div>> <span style="white-space:pre-wrap">  </span>}</div><div>>       else</div><div>> <span style="white-space:pre-wrap">  </span>{</div><div>> <span style="white-space:pre-wrap">   </span>  middle[0] = PETSC_i;</div><div>> <span style="white-space:pre-wrap">     </span>  middle[1] = drive + randoms_loc[rowcount];</div><div>> <span style="white-space:pre-wrap">       </span>  middle[2] = PETSC_i;</div><div>> <span style="white-space:pre-wrap">     </span>  idxmiddle[0] = rowcount - 1;</div><div>> <span style="white-space:pre-wrap">     </span>  idxmiddle[1] = rowcount;</div><div>> <span style="white-space:pre-wrap"> </span>  idxmiddle[2] = rowcount + 1;</div><div>> <span style="white-space:pre-wrap">     </span>  ierr =</div><div>> <span style="white-space:pre-wrap">   </span>    MatSetValues (A, 1, &rowcount, 3, idxmiddle, middle,</div><div>> <span style="white-space:pre-wrap">                      </span>  INSERT_VALUES);</div><div>> <span style="white-space:pre-wrap">  </span>  CHKERRQ (ierr);</div><div>> <span style="white-space:pre-wrap">  </span>}</div></div></div><div>683a625,626</div><span class=""><div>>   ierr = VecRestoreArray (appctx->randoms, &randoms_loc);</div><div>>   CHKERRQ (ierr);</div></span><div>689,691c632,634</div><span class=""><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></span><div>693,704c636,639</div><span class=""><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></span></div><div><br></div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><span><font color="#888888"><div><br></div>-- <br><span class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener
</span></font></span></div></div>
</blockquote></div><br><br clear="all"><span class=""><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>
</span></div></div>
</blockquote></div><br><br clear="all"><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>