Actual source code: ex18.c
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
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: */
14: #include <petscts.h>
16: PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
17: PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
19: int main(int argc, char **argv)
20: {
21: TS ts;
22: Vec x;
23: PetscBool dae = PETSC_TRUE;
25: PetscFunctionBeginUser;
26: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
27: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dae", &dae, NULL));
29: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
30: PetscCall(TSSetIFunction(ts, NULL, IFunction, &dae));
31: PetscCall(TSSetIJacobian(ts, NULL, NULL, IJacobian, &dae));
33: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
34: PetscCall(VecSetSizes(x, 2, PETSC_DECIDE));
35: PetscCall(VecSetFromOptions(x));
36: PetscCall(VecSetUp(x));
37: PetscCall(VecSet(x, 0.5));
38: PetscCall(TSSetSolution(ts, x));
39: PetscCall(VecDestroy(&x));
41: PetscCall(TSSetTimeStep(ts, 1.0 / 16.0));
42: PetscCall(TSSetMaxTime(ts, 1.0));
43: PetscCall(TSSetFromOptions(ts));
44: PetscCall(TSSolve(ts, NULL));
46: PetscCall(TSDestroy(&ts));
47: PetscCall(PetscFinalize());
48: return 0;
49: }
51: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
52: {
53: const PetscScalar *xdot, *x;
54: PetscScalar *f;
55: PetscBool dae = *(PetscBool *)(ctx);
57: PetscFunctionBeginUser;
58: PetscCall(VecGetArrayRead(Xdot, &xdot));
59: PetscCall(VecGetArrayRead(X, &x));
60: PetscCall(VecGetArrayWrite(F, &f));
61: if (dae) {
62: f[0] = x[0] * xdot[0] - x[0] * x[1];
63: f[1] = x[0] - x[1];
64: } else {
65: f[0] = xdot[0] - x[0];
66: f[1] = xdot[1] - x[1];
67: }
68: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
69: PetscCall(VecRestoreArrayRead(X, &x));
70: PetscCall(VecRestoreArrayWrite(F, &f));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat A, Mat B, void *ctx)
75: {
76: const PetscScalar *xdot, *x;
77: PetscBool dae = *(PetscBool *)(ctx);
79: PetscFunctionBeginUser;
80: PetscCall(VecGetArrayRead(Xdot, &xdot));
81: PetscCall(VecGetArrayRead(X, &x));
82: if (dae) {
83: PetscCall(MatSetValue(B, 0, 0, shift * x[0] + xdot[0] - x[1], INSERT_VALUES));
84: PetscCall(MatSetValue(B, 0, 1, -x[0], INSERT_VALUES));
85: PetscCall(MatSetValue(B, 1, 0, 1.0, INSERT_VALUES));
86: PetscCall(MatSetValue(B, 1, 1, -1.0, INSERT_VALUES));
87: } else {
88: PetscCall(MatZeroEntries(B));
89: PetscCall(MatShift(B, shift));
90: }
91: PetscCall(VecRestoreArrayRead(X, &x));
92: PetscCall(VecRestoreArrayRead(Xdot, &xdot));
93: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
94: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
95: if (A != B) {
96: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
97: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
98: }
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: /*TEST
104: testset:
105: args: -ts_view_solution -ts_max_steps 10 -ts_dt 0.1 -ts_view_solution -ts_adapt_type {{none basic}} -ts_exact_final_time matchstep
107: test:
108: output_file: output/ex18_1.out
109: suffix: bdf
110: args: -ts_type bdf
112: test:
113: output_file: output/ex18_1.out
114: suffix: dirk
115: 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}}
117: TEST*/