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, random = PETSC_FALSE;

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

 30:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
 31:   PetscCall(TSSetIFunction(ts, NULL, IFunction, &dae));
 32:   PetscCall(TSSetIJacobian(ts, NULL, NULL, IJacobian, &dae));

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

 45:   PetscCall(TSSetTimeStep(ts, 1.0 / 16.0));
 46:   PetscCall(TSSetMaxTime(ts, 1.0));
 47:   PetscCall(TSSetFromOptions(ts));
 48:   PetscCall(TSSolve(ts, NULL));

 50:   PetscCall(TSDestroy(&ts));
 51:   PetscCall(PetscFinalize());
 52:   return 0;
 53: }

 55: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
 56: {
 57:   const PetscScalar *xdot, *x;
 58:   PetscScalar       *f;
 59:   PetscBool          dae = *(PetscBool *)(ctx);

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

 78: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat A, Mat B, void *ctx)
 79: {
 80:   const PetscScalar *xdot, *x;
 81:   PetscBool          dae = *(PetscBool *)(ctx);

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

106: /*TEST

108:   testset:
109:     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 -snes_error_if_not_converged

111:     test:
112:       output_file: output/ex18_1.out
113:       suffix: bdf
114:       args: -ts_type bdf

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

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

126: TEST*/