Bug Summary

File:ex74.c
Warning:line 239, column 3
Value stored to 'ierr' is never read

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1static char help[] = "\
2Solves the constant-coefficient 1D Heat equation with an Implicit \n\
3Runge-Kutta method using MatKAIJ. \n\
4 \n\
5 du d^2 u \n\
6 -- = a ----- ; 0 <= x <= 1; \n\
7 dt dx^2 \n\
8 \n\
9 with periodic boundary conditions \n\
10 \n\
112nd order central discretization in space: \n\
12 \n\
13 [ d^2 u ] u_{i+1} - 2u_i + u_{i-1} \n\
14 [ ----- ] = ------------------------ \n\
15 [ dx^2 ]i h^2 \n\
16 \n\
17 i = grid index; h = x_{i+1}-x_i (Uniform) \n\
18 0 <= i < n h = 1.0/n \n\
19 \n\
20Thus, \n\
21 \n\
22 du \n\
23 -- = Ju; J = (a/h^2) tridiagonal(1,-2,1)_n \n\
24 dt \n\
25 \n\
26Implicit Runge-Kutta method: \n\
27 \n\
28 U^(k) = u^n + dt \\sum_i a_{ki} JU^{i} \n\
29 u^{n+1} = u^n + dt \\sum_i b_i JU^{i} \n\
30 \n\
31 i = 1,...,s (s -> number of stages) \n\
32 \n\
33At each time step, we solve \n\
34 \n\
35 [ 1 ] 1 \n\
36 [ -- I \\otimes A^{-1} - J \\otimes I ] U = -- u^n \\otimes A^{-1} \n\
37 [ dt ] dt \n\
38 \n\
39 where A is the Butcher tableaux of the implicit \n\
40 Runge-Kutta method, \n\
41 \n\
42with MATKAIJ and KSP. \n\
43 \n\
44Available IRK Methods: \n\
45 gauss n-stage Gauss method \n\
46 \n";
47
48/*T
49 Concepts: MATKAIJ
50 Concepts: MAT
51 Concepts: KSP
52T*/
53
54/*
55 Include "petscksp.h" so that we can use KSP solvers. Note that this file
56 automatically includes:
57 petscsys.h - base PETSc routines
58 petscvec.h - vectors
59 petscmat.h - matrices
60 petscis.h - index sets
61 petscviewer.h - viewers
62 petscpc.h - preconditioners
63*/
64#include <petscksp.h>
65#include <petscdt.h>
66
67/* define the IRK methods available */
68#define IRKGAUSS"gauss" "gauss"
69
70typedef enum {
71 PHYSICS_DIFFUSION,
72 PHYSICS_ADVECTION
73} PhysicsType;
74const char *const PhysicsTypes[] = {"DIFFUSION","ADVECTION","PhysicsType","PHYSICS_",NULL((void*)0)};
75
76typedef struct __context__ {
77 PetscReal a; /* diffusion coefficient */
78 PetscReal xmin,xmax; /* domain bounds */
79 PetscInt imax; /* number of grid points */
80 PetscInt niter; /* number of time iterations */
81 PetscReal dt; /* time step size */
82 PhysicsType physics_type;
83} UserContext;
84
85static PetscErrorCode ExactSolution(Vec,void*,PetscReal);
86static PetscErrorCode RKCreate_Gauss(PetscInt,PetscScalar**,PetscScalar**,PetscReal**);
87static PetscErrorCode Assemble_AdvDiff(MPI_Comm,UserContext*,Mat*);
88
89#include <petsc/private/kernels/blockinvert.h>
90
91int main(int argc, char **argv)
92{
93 PetscErrorCode ierr;
94 Vec u,uex,rhs,z;
95 UserContext ctxt;
96 PetscInt nstages,is,ie,matis,matie,*ix,*ix2;
97 PetscInt n,i,s,t,total_its;
98 PetscScalar *A,*B,*At,*b,*zvals,one = 1.0;
99 PetscReal *c,err,time;
100 Mat Identity,J,TA,SC,R;
101 KSP ksp;
102 PetscFunctionList IRKList = NULL((void*)0);
103 char irktype[256] = IRKGAUSS"gauss";
104
105 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
106 ierr = PetscFunctionListAdd(&IRKList,IRKGAUSS,RKCreate_Gauss)PetscFunctionListAdd_Private((&IRKList),("gauss"),(PetscVoidFunction
)(RKCreate_Gauss));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),106,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
107
108 /* default value */
109 ctxt.a = 1.0;
110 ctxt.xmin = 0.0;
111 ctxt.xmax = 1.0;
112 ctxt.imax = 20;
113 ctxt.niter = 0;
114 ctxt.dt = 0.0;
115 ctxt.physics_type = PHYSICS_DIFFUSION;
116
117 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"IRK options","")0; do { PetscOptionItems PetscOptionsObjectBase; PetscOptionItems
*PetscOptionsObject = &PetscOptionsObjectBase; PetscMemzero
(PetscOptionsObject,sizeof(PetscOptionItems)); for (PetscOptionsObject
->count=(PetscOptionsPublish?-1:1); PetscOptionsObject->
count<2; PetscOptionsObject->count++) { PetscErrorCode _5_ierr
= PetscOptionsBegin_Private(PetscOptionsObject,PETSC_COMM_WORLD
,((void*)0),"IRK options","");do {if (__builtin_expect(!!(_5_ierr
),0)) return PetscError(((MPI_Comm)0x44000001),117,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),117,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
118 ierr = PetscOptionsReal("-a","diffusion coefficient","<1.0>",ctxt.a,&ctxt.a,NULL)PetscOptionsReal_Private(PetscOptionsObject,"-a","diffusion coefficient"
,"<1.0>",ctxt.a,&ctxt.a,((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),118,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
119 ierr = PetscOptionsInt ("-imax","grid size","<20>",ctxt.imax,&ctxt.imax,NULL)PetscOptionsInt_Private(PetscOptionsObject,"-imax","grid size"
,"<20>",ctxt.imax,&ctxt.imax,((void*)0),(-2147483647
- 1),2147483647);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),119,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
120 ierr = PetscOptionsReal("-xmin","xmin","<0.0>",ctxt.xmin,&ctxt.xmin,NULL)PetscOptionsReal_Private(PetscOptionsObject,"-xmin","xmin","<0.0>"
,ctxt.xmin,&ctxt.xmin,((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),120,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
121 ierr = PetscOptionsReal("-xmax","xmax","<1.0>",ctxt.xmax,&ctxt.xmax,NULL)PetscOptionsReal_Private(PetscOptionsObject,"-xmax","xmax","<1.0>"
,ctxt.xmax,&ctxt.xmax,((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),121,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
122 ierr = PetscOptionsInt ("-niter","number of time steps","<0>",ctxt.niter,&ctxt.niter,NULL)PetscOptionsInt_Private(PetscOptionsObject,"-niter","number of time steps"
,"<0>",ctxt.niter,&ctxt.niter,((void*)0),(-2147483647
- 1),2147483647);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),122,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
123 ierr = PetscOptionsReal("-dt","time step size","<0.0>",ctxt.dt,&ctxt.dt,NULL)PetscOptionsReal_Private(PetscOptionsObject,"-dt","time step size"
,"<0.0>",ctxt.dt,&ctxt.dt,((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),123,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
124 ierr = PetscOptionsFList("-irk_type","IRK method family","",IRKList,irktype,irktype,sizeof(irktype),NULL)PetscOptionsFList_Private(PetscOptionsObject,"-irk_type","IRK method family"
,"",IRKList,irktype,irktype,sizeof(irktype),((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),124,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
125 nstages = 2;
126 ierr = PetscOptionsInt ("-irk_nstages","Number of stages in IRK method","",nstages,&nstages,NULL)PetscOptionsInt_Private(PetscOptionsObject,"-irk_nstages","Number of stages in IRK method"
,"",nstages,&nstages,((void*)0),(-2147483647 - 1),2147483647
);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),126,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
127 ierr = PetscOptionsEnum("-physics_type","Type of process to discretize","",PhysicsTypes,(PetscEnum)ctxt.physics_type,(PetscEnum*)&ctxt.physics_type,NULL)PetscOptionsEnum_Private(PetscOptionsObject,"-physics_type","Type of process to discretize"
,"",PhysicsTypes,(PetscEnum)ctxt.physics_type,(PetscEnum*)&
ctxt.physics_type,((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),127,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
128 ierr = PetscOptionsEnd()_5_ierr = PetscOptionsEnd_Private(PetscOptionsObject);do {if (
__builtin_expect(!!(_5_ierr),0)) return PetscError(((MPI_Comm
)0x44000001),128,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,_5_ierr,PETSC_ERROR_REPEAT," ");} while (0);}} while (0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),128,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
129
130 /* allocate and initialize solution vector and exact solution */
131 ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),131,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
132 ierr = VecSetSizes(u,PETSC_DECIDE-1,ctxt.imax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),132,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
133 ierr = VecSetFromOptions(u);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),133,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
134 ierr = VecDuplicate(u,&uex);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),134,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
135 /* initial solution */
136 ierr = ExactSolution(u ,&ctxt,0.0);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),136,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
137 /* exact solution */
138 ierr = ExactSolution(uex,&ctxt,ctxt.dt*ctxt.niter);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),138,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
139
140 { /* Create A,b,c */
141 PetscErrorCode (*irkcreate)(PetscInt,PetscScalar**,PetscScalar**,PetscReal**);
142 ierr = PetscFunctionListFind(IRKList,irktype,&irkcreate)PetscFunctionListFind_Private((IRKList),(irktype),(PetscVoidFunction
*)(&irkcreate));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),142,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
143 ierr = (*irkcreate)(nstages,&A,&b,&c);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),143,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
144 }
145 { /* Invert A */
146 /* PETSc does not provide a routine to calculate the inverse of a general matrix.
147 * To get the inverse of A, we form a sequential BAIJ matrix from it, consisting of a single block with block size
148 * equal to the dimension of A, and then use MatInvertBlockDiagonal(). */
149 Mat A_baij;
150 PetscInt idxm[1]={0},idxn[1]={0};
151 const PetscScalar *A_inv;
152 ierr = MatCreateSeqBAIJ(PETSC_COMM_SELF((MPI_Comm)0x44000001),nstages,nstages,nstages,1,NULL((void*)0),&A_baij);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),152,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
153 ierr = MatSetOption(A_baij,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),153,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
154 ierr = MatSetValuesBlocked(A_baij,1,idxm,1,idxn,A,INSERT_VALUES);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),154,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
155 ierr = MatAssemblyBegin(A_baij,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),155,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
156 ierr = MatAssemblyEnd(A_baij,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),156,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
157 ierr = MatInvertBlockDiagonal(A_baij,&A_inv);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),157,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
158 ierr = PetscMemcpy(A,A_inv,nstages*nstages*sizeof(PetscScalar));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),158,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
159 ierr = MatDestroy(&A_baij);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),159,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
160 }
161 /* Scale (1/dt)*A^{-1} and (1/dt)*b */
162 for (s=0; s<nstages*nstages; s++) A[s] *= 1.0/ctxt.dt;
163 for (s=0; s<nstages; s++) b[s] *= (-ctxt.dt);
164
165 /* Compute row sums At and identity B */
166 ierr = PetscMalloc2(nstages,&At,PetscSqr(nstages),&B)PetscMallocA(2,PETSC_FALSE,166,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,(size_t)(nstages)*sizeof(**(&At)),(&At),(size_t)(((nstages
)*(nstages)))*sizeof(**(&B)),(&B));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),166,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
167 for (s=0; s<nstages; s++) {
168 At[s] = 0;
169 for (t=0; t<nstages; t++) {
170 At[s] += A[s+nstages*t]; /* Row sums of */
171 B[s+nstages*t] = 1.*(s == t); /* identity */
172 }
173 }
174
175 /* allocate and calculate the (-J) matrix */
176 switch (ctxt.physics_type) {
177 case PHYSICS_ADVECTION:
178 case PHYSICS_DIFFUSION:
179 ierr = Assemble_AdvDiff(PETSC_COMM_WORLD,&ctxt,&J);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),179,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
180 }
181 ierr = MatCreate(PETSC_COMM_WORLD,&Identity);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),181,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
182 ierr = MatSetType(Identity,MATAIJ"aij");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),182,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
183 ierr = MatGetOwnershipRange(J,&matis,&matie);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),183,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
184 ierr = MatSetSizes(Identity,matie-matis,matie-matis,ctxt.imax,ctxt.imax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),184,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
185 ierr = MatSetUp(Identity);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),185,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
186 for (i=matis; i<matie; i++) {
187 ierr= MatSetValues(Identity,1,&i,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),187,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
188 }
189 ierr = MatAssemblyBegin(Identity,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),189,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
190 ierr = MatAssemblyEnd (Identity,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),190,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
191
192 /* Create the KAIJ matrix for solving the stages */
193 ierr = MatCreateKAIJ(J,nstages,nstages,A,B,&TA);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),193,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
194
195 /* Create the KAIJ matrix for step completion */
196 ierr = MatCreateKAIJ(J,1,nstages,NULL((void*)0),b,&SC);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),196,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
197
198 /* Create the KAIJ matrix to create the R for solving the stages */
199 ierr = MatCreateKAIJ(Identity,nstages,1,NULL((void*)0),At,&R);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),199,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
200
201 /* Create and set options for KSP */
202 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),202,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
203 ierr = KSPSetOperators(ksp,TA,TA);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),203,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
204 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),204,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
205
206 /* Allocate work and right-hand-side vectors */
207 ierr = VecCreate(PETSC_COMM_WORLD,&z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),207,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
208 ierr = VecSetFromOptions(z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),208,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
209 ierr = VecSetSizes(z,PETSC_DECIDE-1,ctxt.imax*nstages);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),209,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
210 ierr = VecDuplicate(z,&rhs);
211
212 ierr = VecGetOwnershipRange(u,&is,&ie);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),212,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
213 ierr = PetscMalloc3(nstages,&ix,nstages,&zvals,ie-is,&ix2)PetscMallocA(3,PETSC_FALSE,213,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,(size_t)(nstages)*sizeof(**(&ix)),(&ix),(size_t)(nstages
)*sizeof(**(&zvals)),(&zvals),(size_t)(ie-is)*sizeof(
**(&ix2)),(&ix2));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),213,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
214 /* iterate in time */
215 for (n=0,time=0.,total_its=0; n<ctxt.niter; n++) {
216 PetscInt its;
217
218 /* compute and set the right hand side */
219 ierr = MatMult(R,u,rhs);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),219,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
220
221 /* Solve the system */
222 ierr = KSPSolve(ksp,rhs,z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),222,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
223 ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),223,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
224 total_its += its;
225
226 /* Update the solution */
227 ierr = MatMultAdd(SC,z,u,u);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),227,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
228
229 /* time step complete */
230 time += ctxt.dt;
231 }
232 PetscFree3(ix,ix2,zvals)PetscFreeA(3,232,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,&(ix),&(ix2),&(zvals));
233
234 /* Deallocate work and right-hand-side vectors */
235 ierr = VecDestroy(&z);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),235,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
236 ierr = VecDestroy(&rhs);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),236,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
237
238 /* Calculate error in final solution */
239 ierr = VecAYPX(uex,-1.0,u);
Value stored to 'ierr' is never read
240 ierr = VecNorm(uex,NORM_2,&err);
241 err = PetscSqrtReal(err*err/((PetscReal)ctxt.imax))sqrt(err*err/((PetscReal)ctxt.imax));
242 ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 norm of the numerical error = %g (time=%g)\n",(double)err,(double)time);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),242,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
243 ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of time steps: %D (%D Krylov iterations)\n",ctxt.niter,total_its);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),243,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
244
245 /* Free up memory */
246 ierr = KSPDestroy(&ksp); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),246,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
247 ierr = MatDestroy(&TA); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),247,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
248 ierr = MatDestroy(&SC); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),248,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
249 ierr = MatDestroy(&R); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),249,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
250 ierr = MatDestroy(&J); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),250,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
251 ierr = MatDestroy(&Identity); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),251,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
252 ierr = PetscFree3(A,b,c)PetscFreeA(3,252,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,&(A),&(b),&(c)); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),252,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
253 ierr = PetscFree2(At,B)PetscFreeA(2,253,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,&(At),&(B)); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),253,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
254 ierr = VecDestroy(&uex); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),254,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
255 ierr = VecDestroy(&u); CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),255,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
256 ierr = PetscFunctionListDestroy(&IRKList);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),256,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
257
258 ierr = PetscFinalize();
259 return ierr;
260}
261
262PetscErrorCode ExactSolution(Vec u,void *c,PetscReal t)
263{
264 UserContext *ctxt = (UserContext*) c;
265 PetscErrorCode ierr;
266 PetscInt i,is,ie;
267 PetscScalar *uarr;
268 PetscReal x,dx,a=ctxt->a,pi=PETSC_PI3.1415926535897932384626433832795029;
269
270 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
; petscstack->line[petscstack->currentsize] = 270; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);
271 dx = (ctxt->xmax - ctxt->xmin)/((PetscReal) ctxt->imax);
272 ierr = VecGetOwnershipRange(u,&is,&ie);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),272,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
273 ierr = VecGetArray(u,&uarr);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),273,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
274 for(i=is; i<ie; i++) {
275 x = i * dx;
276 switch (ctxt->physics_type) {
277 case PHYSICS_DIFFUSION:
278 uarr[i-is] = PetscExpScalar(-4.0*pi*pi*a*t)exp(-4.0*pi*pi*a*t)*PetscSinScalar(2*pi*x)sin(2*pi*x);
279 break;
280 case PHYSICS_ADVECTION:
281 uarr[i-is] = PetscSinScalar(2*pi*(x - a*t))sin(2*pi*(x - a*t));
282 break;
283 default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[ctxt->physics_type])return PetscError(((MPI_Comm)0x44000001),283,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,56,PETSC_ERROR_INITIAL,"No support for physics type %s",PhysicsTypes
[ctxt->physics_type]);
284 }
285 }
286 ierr = VecRestoreArray(u,&uarr);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),286,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
287 PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0);
288}
289
290/* Arrays should be freed with PetscFree3(A,b,c) */
291static PetscErrorCode RKCreate_Gauss(PetscInt nstages,PetscScalar **gauss_A,PetscScalar **gauss_b,PetscReal **gauss_c)
292{
293 PetscErrorCode ierr;
294 PetscScalar *A,*G0,*G1;
295 PetscReal *b,*c;
296 PetscInt i,j;
297 Mat G0mat,G1mat,Amat;
298
299 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
; petscstack->line[petscstack->currentsize] = 299; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);
300 ierr = PetscMalloc3(PetscSqr(nstages),&A,nstages,gauss_b,nstages,&c)PetscMallocA(3,PETSC_FALSE,300,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,(size_t)(((nstages)*(nstages)))*sizeof(**(&A)),(&A),
(size_t)(nstages)*sizeof(**(gauss_b)),(gauss_b),(size_t)(nstages
)*sizeof(**(&c)),(&c));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),300,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
301 ierr = PetscMalloc3(nstages,&b,PetscSqr(nstages),&G0,PetscSqr(nstages),&G1)PetscMallocA(3,PETSC_FALSE,301,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,(size_t)(nstages)*sizeof(**(&b)),(&b),(size_t)(((nstages
)*(nstages)))*sizeof(**(&G0)),(&G0),(size_t)(((nstages
)*(nstages)))*sizeof(**(&G1)),(&G1));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),301,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
302 ierr = PetscDTGaussQuadrature(nstages,0.,1.,c,b);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),302,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
303 for (i=0; i<nstages; i++) (*gauss_b)[i] = b[i]; /* copy to possibly-complex array */
304
305 /* A^T = G0^{-1} G1 */
306 for (i=0; i<nstages; i++) {
307 for (j=0; j<nstages; j++) {
308 G0[i*nstages+j] = PetscPowRealInt(c[i],j);
309 G1[i*nstages+j] = PetscPowRealInt(c[i],j+1)/(j+1);
310 }
311 }
312 /* The arrays above are row-aligned, but we create dense matrices as the transpose */
313 ierr = MatCreateSeqDense(PETSC_COMM_SELF((MPI_Comm)0x44000001),nstages,nstages,G0,&G0mat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),313,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
314 ierr = MatCreateSeqDense(PETSC_COMM_SELF((MPI_Comm)0x44000001),nstages,nstages,G1,&G1mat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),314,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
315 ierr = MatCreateSeqDense(PETSC_COMM_SELF((MPI_Comm)0x44000001),nstages,nstages,A,&Amat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),315,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
316 ierr = MatLUFactor(G0mat,NULL((void*)0),NULL((void*)0),NULL((void*)0));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),316,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
317 ierr = MatMatSolve(G0mat,G1mat,Amat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),317,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
318 ierr = MatTranspose(Amat,MAT_INPLACE_MATRIX,&Amat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),318,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
319
320 ierr = MatDestroy(&G0mat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),320,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
321 ierr = MatDestroy(&G1mat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),321,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
322 ierr = MatDestroy(&Amat);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),322,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
323 ierr = PetscFree3(b,G0,G1)PetscFreeA(3,323,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,&(b),&(G0),&(G1));CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),323,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
324 *gauss_A = A;
325 *gauss_c = c;
326 PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0);
327}
328
329static PetscErrorCode Assemble_AdvDiff(MPI_Comm comm,UserContext *user,Mat *J)
330{
331 PetscErrorCode ierr;
332 PetscInt matis,matie,i;
333 PetscReal dx,dx2;
334
335 PetscFunctionBegindo { do { ; if (petscstack && (petscstack->currentsize
< 64)) { petscstack->function[petscstack->currentsize
] = __func__; petscstack->file[petscstack->currentsize]
= "/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
; petscstack->line[petscstack->currentsize] = 335; petscstack
->petscroutine[petscstack->currentsize] = PETSC_TRUE; petscstack
->currentsize++; } if (petscstack) { petscstack->hotdepth
+= (PETSC_FALSE || petscstack->hotdepth); } ; } while (0)
; ; } while (0);
336 dx = (user->xmax - user->xmin)/((PetscReal)user->imax); dx2 = dx*dx;
337 ierr = MatCreate(comm,J);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),337,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
338 ierr = MatSetType(*J,MATAIJ"aij");CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),338,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
339 ierr = MatSetSizes(*J,PETSC_DECIDE-1,PETSC_DECIDE-1,user->imax,user->imax);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),339,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
340 ierr = MatSetUp(*J);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),340,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
341 ierr = MatGetOwnershipRange(*J,&matis,&matie);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),341,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
342 for (i=matis; i<matie; i++) {
343 PetscScalar values[3];
344 PetscInt col[3];
345 switch (user->physics_type) {
346 case PHYSICS_DIFFUSION:
347 values[0] = -user->a*1.0/dx2;
348 values[1] = user->a*2.0/dx2;
349 values[2] = -user->a*1.0/dx2;
350 break;
351 case PHYSICS_ADVECTION:
352 values[0] = -user->a*.5/dx;
353 values[1] = 0.;
354 values[2] = user->a*.5/dx;
355 break;
356 default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for physics type %s",PhysicsTypes[user->physics_type])return PetscError(((MPI_Comm)0x44000001),356,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,56,PETSC_ERROR_INITIAL,"No support for physics type %s",PhysicsTypes
[user->physics_type]);
357 }
358 /* periodic boundaries */
359 if (i == 0) {
360 col[0] = user->imax-1;
361 col[1] = i;
362 col[2] = i+1;
363 } else if (i == user->imax-1) {
364 col[0] = i-1;
365 col[1] = i;
366 col[2] = 0;
367 } else {
368 col[0] = i-1;
369 col[1] = i;
370 col[2] = i+1;
371 }
372 ierr= MatSetValues(*J,1,&i,3,col,values,INSERT_VALUES);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),372,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
373 }
374 ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),374,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
375 ierr = MatAssemblyEnd (*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr)do {if (__builtin_expect(!!(ierr),0)) return PetscError(((MPI_Comm
)0x44000001),375,__func__,"/sandbox/petsc/petsc.next/src/ksp/ksp/examples/tutorials/ex74.c"
,ierr,PETSC_ERROR_REPEAT," ");} while (0);
376 PetscFunctionReturn(0)do { do { ; if (petscstack && petscstack->currentsize
> 0) { petscstack->currentsize--; petscstack->function
[petscstack->currentsize] = 0; petscstack->file[petscstack
->currentsize] = 0; petscstack->line[petscstack->currentsize
] = 0; petscstack->petscroutine[petscstack->currentsize
] = PETSC_FALSE; } if (petscstack) { petscstack->hotdepth =
(((petscstack->hotdepth-1)<(0)) ? (0) : (petscstack->
hotdepth-1)); } ; } while (0); return(0);} while (0);
377}
378
379/*TEST
380 test:
381 suffix: 1
382 args: -a 0.1 -dt .125 -niter 5 -imax 40 -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -irk_type gauss -irk_nstages 2
383 test:
384 suffix: 2
385 args: -a 0.1 -dt .125 -niter 5 -imax 40 -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -irk_type gauss -irk_nstages 4 -ksp_gmres_restart 100
386 test:
387 suffix: 3
388 args: -a 1 -dt .33 -niter 3 -imax 40 -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -irk_type gauss -irk_nstages 4 -ksp_gmres_restart 100 -physics_type advection
389
390TEST*/