[petsc-users] TS example tutorial example 4 adapted - not parallelizing
Matthew Knepley
knepley at gmail.com
Wed Oct 15 05:45:28 CDT 2014
On Wed, Oct 15, 2014 at 3:20 AM, Analabha Roy <hariseldon99 at gmail.com>
wrote:
>
> 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?
>
I am not sure what you are doing at line 483. You seem to have all
processes indexing a Vec
at all entries, which cannot be correct.
Matt
> 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/
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141015/fc9b0dac/attachment-0001.html>
More information about the petsc-users
mailing list