1 static char help[] = "Solves the constant-coefficient 1D heat equation \n\
2 with an Implicit Runge-Kutta method using MatKAIJ. \n\
3 \n\
4 du d^2 u \n\
5 -- = a ----- ; 0 <= x <= 1; \n\
6 dt dx^2 \n\
7 \n\
8 with periodic boundary conditions \n\
9 \n\
10 2nd order central discretization in space: \n\
11 \n\
12 [ d^2 u ] u_{i+1} - 2u_i + u_{i-1} \n\
13 [ ----- ] = ------------------------ \n\
14 [ dx^2 ]i h^2 \n\
15 \n\
16 i = grid index; h = x_{i+1}-x_i (Uniform) \n\
17 0 <= i < n h = 1.0/n \n\
18 \n\
19 Thus, \n\
20 \n\
21 du \n\
22 -- = Ju; J = (a/h^2) tridiagonal(1,-2,1)_n \n\
23 dt \n\
24 \n\
25 This example is a TS version of the KSP ex74.c tutorial. \n";
26
27 #include <petscts.h>
28
29 typedef enum {
30 PHYSICS_DIFFUSION,
31 PHYSICS_ADVECTION
32 } PhysicsType;
33 const char *const PhysicsTypes[] = {"DIFFUSION", "ADVECTION", "PhysicsType", "PHYSICS_", NULL};
34
35 typedef struct Context {
36 PetscReal a; /* diffusion coefficient */
37 PetscReal xmin, xmax; /* domain bounds */
38 PetscInt imax; /* number of grid points */
39 PhysicsType physics_type;
40 } UserContext;
41
42 static PetscErrorCode ExactSolution(Vec, void *, PetscReal);
43 static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
44
main(int argc,char ** argv)45 int main(int argc, char **argv)
46 {
47 TS ts;
48 Mat A;
49 Vec u, uex;
50 UserContext ctxt;
51 PetscReal err, ftime;
52
53 PetscFunctionBeginUser;
54 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
55 /* default value */
56 ctxt.a = 0.1;
57 ctxt.xmin = 0.0;
58 ctxt.xmax = 1.0;
59 ctxt.imax = 40;
60 ctxt.physics_type = PHYSICS_DIFFUSION;
61
62 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "IRK options", "");
63 PetscCall(PetscOptionsReal("-a", "diffusion coefficient", "<1.0>", ctxt.a, &ctxt.a, NULL));
64 PetscCall(PetscOptionsInt("-imax", "grid size", "<20>", ctxt.imax, &ctxt.imax, NULL));
65 PetscCall(PetscOptionsReal("-xmin", "xmin", "<0.0>", ctxt.xmin, &ctxt.xmin, NULL));
66 PetscCall(PetscOptionsReal("-xmax", "xmax", "<1.0>", ctxt.xmax, &ctxt.xmax, NULL));
67 PetscCall(PetscOptionsEnum("-physics_type", "Type of process to discretize", "", PhysicsTypes, (PetscEnum)ctxt.physics_type, (PetscEnum *)&ctxt.physics_type, NULL));
68 PetscOptionsEnd();
69
70 /* allocate and initialize solution vector and exact solution */
71 PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
72 PetscCall(VecSetSizes(u, PETSC_DECIDE, ctxt.imax));
73 PetscCall(VecSetFromOptions(u));
74 PetscCall(VecDuplicate(u, &uex));
75 /* initial solution */
76 PetscCall(ExactSolution(u, &ctxt, 0.0));
77
78 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
79 PetscCall(MatSetType(A, MATAIJ));
80 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, ctxt.imax, ctxt.imax));
81 PetscCall(MatSetUp(A));
82
83 /* Create and set options for TS */
84 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
85 PetscCall(TSSetProblemType(ts, TS_LINEAR));
86 PetscCall(TSSetTimeStep(ts, 0.125));
87 PetscCall(TSSetSolution(ts, u));
88 PetscCall(TSSetMaxSteps(ts, 10));
89 PetscCall(TSSetMaxTime(ts, 1.0));
90 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
91 PetscCall(TSSetRHSFunction(ts, NULL, TSComputeRHSFunctionLinear, &ctxt));
92 PetscCall(RHSJacobian(ts, 0, u, A, A, &ctxt));
93 PetscCall(TSSetRHSJacobian(ts, A, A, TSComputeRHSJacobianConstant, &ctxt));
94 PetscCall(TSSetFromOptions(ts));
95 PetscCall(TSSolve(ts, u));
96
97 PetscCall(TSGetSolveTime(ts, &ftime));
98 /* exact solution */
99 PetscCall(ExactSolution(uex, &ctxt, ftime));
100
101 /* Calculate error in final solution */
102 PetscCall(VecAYPX(uex, -1.0, u));
103 PetscCall(VecNorm(uex, NORM_2, &err));
104 err = PetscSqrtReal(err * err / ((PetscReal)ctxt.imax));
105 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L2 norm of the numerical error = %g (time=%g)\n", (double)err, (double)ftime));
106
107 /* Free up memory */
108 PetscCall(TSDestroy(&ts));
109 PetscCall(MatDestroy(&A));
110 PetscCall(VecDestroy(&uex));
111 PetscCall(VecDestroy(&u));
112 PetscCall(PetscFinalize());
113 return 0;
114 }
115
ExactSolution(Vec u,void * c,PetscReal t)116 PetscErrorCode ExactSolution(Vec u, void *c, PetscReal t)
117 {
118 UserContext *ctxt = (UserContext *)c;
119 PetscInt i, is, ie;
120 PetscScalar *uarr;
121 PetscReal x, dx, a = ctxt->a, pi = PETSC_PI;
122
123 PetscFunctionBeginUser;
124 dx = (ctxt->xmax - ctxt->xmin) / ((PetscReal)ctxt->imax);
125 PetscCall(VecGetOwnershipRange(u, &is, &ie));
126 PetscCall(VecGetArray(u, &uarr));
127 for (i = is; i < ie; i++) {
128 x = i * dx;
129 switch (ctxt->physics_type) {
130 case PHYSICS_DIFFUSION:
131 uarr[i - is] = PetscExpScalar(-4.0 * pi * pi * a * t) * PetscSinScalar(2 * pi * x);
132 break;
133 case PHYSICS_ADVECTION:
134 uarr[i - is] = PetscSinScalar(2 * pi * (x - a * t));
135 break;
136 default:
137 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for physics type %s", PhysicsTypes[ctxt->physics_type]);
138 }
139 }
140 PetscCall(VecRestoreArray(u, &uarr));
141 PetscFunctionReturn(PETSC_SUCCESS);
142 }
143
RHSJacobian(TS ts,PetscReal t,Vec U,Mat J,Mat Jpre,PetscCtx ctx)144 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, PetscCtx ctx)
145 {
146 UserContext *user = (UserContext *)ctx;
147 PetscInt matis, matie, i;
148 PetscReal dx, dx2;
149
150 PetscFunctionBeginUser;
151 dx = (user->xmax - user->xmin) / ((PetscReal)user->imax);
152 dx2 = dx * dx;
153 PetscCall(MatGetOwnershipRange(J, &matis, &matie));
154 for (i = matis; i < matie; i++) {
155 PetscScalar values[3];
156 PetscInt col[3];
157 switch (user->physics_type) {
158 case PHYSICS_DIFFUSION:
159 values[0] = user->a * 1.0 / dx2;
160 values[1] = -user->a * 2.0 / dx2;
161 values[2] = user->a * 1.0 / dx2;
162 break;
163 case PHYSICS_ADVECTION:
164 values[0] = user->a * .5 / dx;
165 values[1] = 0.;
166 values[2] = -user->a * .5 / dx;
167 break;
168 default:
169 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for physics type %s", PhysicsTypes[user->physics_type]);
170 }
171 /* periodic boundaries */
172 if (i == 0) {
173 col[0] = user->imax - 1;
174 col[1] = i;
175 col[2] = i + 1;
176 } else if (i == user->imax - 1) {
177 col[0] = i - 1;
178 col[1] = i;
179 col[2] = 0;
180 } else {
181 col[0] = i - 1;
182 col[1] = i;
183 col[2] = i + 1;
184 }
185 PetscCall(MatSetValues(J, 1, &i, 3, col, values, INSERT_VALUES));
186 }
187 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
188 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
189 PetscFunctionReturn(PETSC_SUCCESS);
190 }
191
192 /*TEST
193
194 test:
195 requires: double
196 suffix: 1
197 nsize: {{1 2}}
198 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 2
199
200 test:
201 requires: double
202 suffix: 2
203 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-6 -ts_type irk -ts_irk_nstages 3
204
205 testset:
206 requires: hpddm
207 args: -ts_max_steps 5 -ts_monitor -ksp_monitor_short -pc_type pbjacobi -ksp_atol 1e-4 -ts_type irk -ts_irk_nstages 3 -ksp_view_final_residual -ksp_hpddm_type gcrodr -ksp_type hpddm
208 test:
209 suffix: 3
210 requires: double
211 args: -ksp_hpddm_precision {{single double}shared output}
212 test:
213 suffix: 3_single
214 requires: single
215 args: -ksp_hpddm_precision {{single double}shared output}
216 test:
217 suffix: 3___float128
218 requires: __float128
219 output_file: output/ex74_3.out
220 args: -ksp_hpddm_precision {{double __float128}shared output}
221 TEST*/
222