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 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 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 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