<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Oct 16, 2014 at 4:08 AM, 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> 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></blockquote><div><br></div><div>I am glad its working. One way to find out of bounds access is to use valgrind, which works great.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>I have fixed it. Basically, line 407 changes from the local index (i-ibase) to the global index (just i).</div><span class=""><div><br></div><div> for (i = mybase; i < myend; i++)</div></span><div> if (i == floor ((appctx->m) / 2))</div><span class=""><div> u_localptr[i - mybase] = one;</div><div> else</div><div> u_localptr[i - mybase] = zero; </div><div><br></div></span><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="HOEnZb"><div class="h5"><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>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><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><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><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><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><div>< Compute the exact solution</div><div>---</div></span><span><div>> Compute and store norm of wavefunction in norm_2 </div></span><div>520c460</div><span><div>< ierr = ExactSolution (time, appctx->solution, appctx);</div><div>---</div></span><span><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><div>< /*</div><div>< Print debugging information if desired</div><div>< */</div><div>< if (appctx->debug)</div></span><span><div>< {</div><div>< ierr = PetscPrintf (appctx->comm, "Computed solution vector\n");</div></span><span><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><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><div>> /* Build the lattice sites */</div></span><div>> Vec uabs;</div><span><div>> ierr = VecDuplicate (u, &uabs);</div></span><div>543c469,470</div><span><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><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><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><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><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><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><div>> ierr = VecDestroy (&uabs);</div><div>> CHKERRQ (ierr);</div></span><div>583c508</div><span><div>< #define __FUNCT__ "RHSMatrixHeat"</div><div>---</div><div>> #define __FUNCT__ "RHSMatrixSchrodinger"</div></span><div>585c510</div><span><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><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><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><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><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><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><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><div>> ierr = VecRestoreArray (appctx->randoms, &randoms_loc);</div><div>> CHKERRQ (ierr);</div></span><div>689,691c632,634</div><span><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><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>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><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>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br>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
</div></div>