Actual source code: ex74.c
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";
27: #include <petscts.h>
29: typedef enum {
30: PHYSICS_DIFFUSION,
31: PHYSICS_ADVECTION
32: } PhysicsType;
33: const char *const PhysicsTypes[] = {"DIFFUSION", "ADVECTION", "PhysicsType", "PHYSICS_", NULL};
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;
42: static PetscErrorCode ExactSolution(Vec, void *, PetscReal);
43: static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
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;
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;
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();
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));
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));
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));
97: PetscCall(TSGetSolveTime(ts, &ftime));
98: /* exact solution */
99: PetscCall(ExactSolution(uex, &ctxt, ftime));
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));
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: }
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;
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: }
144: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, void *ctx)
145: {
146: UserContext *user = (UserContext *)ctx;
147: PetscInt matis, matie, i;
148: PetscReal dx, dx2;
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: }
192: /*TEST
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
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
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 quadruple}shared output}
221: TEST*/