1 static char help[] = "Solves a DAE with a non-trivial mass matrix. \n\n";
2 /*
3 Solves:
4 U * dU/dt = U*V
5 U - V = 0
6
7 that can be rewritten in implicit form F = 0, with F as
8 x[0] * xdot[0] - x[0] * x[1]
9 F(t,x,xdot) =
10 x[0] - x[1]
11 It is equivalent to solve dU/dt = U, U = U0 with solution U = U0 * exp(tfinal)
12 */
13
14 #include <petscts.h>
15
16 PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
17 PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
18
main(int argc,char ** argv)19 int main(int argc, char **argv)
20 {
21 TS ts;
22 Vec x;
23 PetscBool dae = PETSC_TRUE, random = PETSC_FALSE;
24
25 PetscFunctionBeginUser;
26 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
27 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dae", &dae, NULL));
28 PetscCall(PetscOptionsGetBool(NULL, NULL, "-random", &random, NULL));
29
30 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
31 PetscCall(TSSetIFunction(ts, NULL, IFunction, &dae));
32 PetscCall(TSSetIJacobian(ts, NULL, NULL, IJacobian, &dae));
33
34 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
35 PetscCall(VecSetSizes(x, 2, PETSC_DECIDE));
36 PetscCall(VecSetFromOptions(x));
37 PetscCall(VecSetUp(x));
38 if (random) {
39 PetscCall(VecSetRandom(x, NULL));
40 PetscCall(VecRealPart(x));
41 } else PetscCall(VecSet(x, 0.5));
42 PetscCall(TSSetSolution(ts, x));
43 PetscCall(VecDestroy(&x));
44
45 PetscCall(TSSetTimeStep(ts, 1.0 / 16.0));
46 PetscCall(TSSetMaxTime(ts, 1.0));
47 PetscCall(TSSetFromOptions(ts));
48 PetscCall(TSSolve(ts, NULL));
49
50 PetscCall(TSDestroy(&ts));
51 PetscCall(PetscFinalize());
52 return 0;
53 }
54
IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)55 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
56 {
57 const PetscScalar *xdot, *x;
58 PetscScalar *f;
59 PetscBool dae = *(PetscBool *)(ctx);
60
61 PetscFunctionBeginUser;
62 PetscCall(VecGetArrayRead(Xdot, &xdot));
63 PetscCall(VecGetArrayRead(X, &x));
64 PetscCall(VecGetArrayWrite(F, &f));
65 if (dae) {
66 f[0] = x[0] * xdot[0] - x[0] * x[1];
67 f[1] = x[0] - x[1];
68 } else {
69 f[0] = xdot[0] - x[0];
70 f[1] = xdot[1] - x[1];
71 }
72 PetscCall(VecRestoreArrayRead(Xdot, &xdot));
73 PetscCall(VecRestoreArrayRead(X, &x));
74 PetscCall(VecRestoreArrayWrite(F, &f));
75 PetscFunctionReturn(PETSC_SUCCESS);
76 }
77
IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal shift,Mat A,Mat B,PetscCtx ctx)78 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat A, Mat B, PetscCtx ctx)
79 {
80 const PetscScalar *xdot, *x;
81 PetscBool dae = *(PetscBool *)(ctx);
82
83 PetscFunctionBeginUser;
84 PetscCall(VecGetArrayRead(Xdot, &xdot));
85 PetscCall(VecGetArrayRead(X, &x));
86 if (dae) {
87 PetscCall(MatSetValue(B, 0, 0, shift * x[0] + xdot[0] - x[1], INSERT_VALUES));
88 PetscCall(MatSetValue(B, 0, 1, -x[0], INSERT_VALUES));
89 PetscCall(MatSetValue(B, 1, 0, 1.0, INSERT_VALUES));
90 PetscCall(MatSetValue(B, 1, 1, -1.0, INSERT_VALUES));
91 } else {
92 PetscCall(MatZeroEntries(B));
93 PetscCall(MatShift(B, shift));
94 }
95 PetscCall(VecRestoreArrayRead(X, &x));
96 PetscCall(VecRestoreArrayRead(Xdot, &xdot));
97 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
98 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
99 if (A != B) {
100 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
101 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
102 }
103 PetscFunctionReturn(PETSC_SUCCESS);
104 }
105
106 /*TEST
107
108 testset:
109 args: -ts_view_solution -ts_max_steps 10 -ts_time_step 0.1 -ts_view_solution -ts_adapt_type {{none basic}} -ts_exact_final_time matchstep -snes_error_if_not_converged
110
111 test:
112 output_file: output/ex18_1.out
113 suffix: bdf
114 args: -ts_type bdf
115
116 test:
117 output_file: output/ex18_1.out
118 suffix: dirk
119 args: -dae {{0 1}} -ts_type dirk -ts_dirk_type {{s212 es122sal es213sal es324sal es325sal 657a es648sa 658a s659a 7510sal es7510sa 759a s7511sal 8614a 8616sal es8516sal}}
120
121 test:
122 output_file: output/ex18_1.out
123 suffix: dirk_explicit_first_random_dae
124 args: -dae -ts_type dirk -ts_dirk_type es122sal -random 1 -ts_max_step_rejections -1
125
126 TEST*/
127