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