1 static char help[] = "Solves the trival ODE 2 du/dt = 1, u(0) = 0. \n\n"; 2 3 #include <petscts.h> 4 #include <petscpc.h> 5 6 PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*); 7 PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); 8 9 int main(int argc,char **argv) 10 { 11 TS ts; 12 Vec x; 13 Vec f; 14 Mat A; 15 PetscErrorCode ierr; 16 17 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 18 19 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 20 ierr = TSSetEquationType(ts,TS_EQ_ODE_IMPLICIT);CHKERRQ(ierr); 21 ierr = VecCreate(PETSC_COMM_WORLD,&f);CHKERRQ(ierr); 22 ierr = VecSetSizes(f,1,PETSC_DECIDE);CHKERRQ(ierr); 23 ierr = VecSetFromOptions(f);CHKERRQ(ierr); 24 ierr = VecSetUp(f);CHKERRQ(ierr); 25 ierr = TSSetIFunction(ts,f,IFunction,NULL);CHKERRQ(ierr); 26 ierr = VecDestroy(&f);CHKERRQ(ierr); 27 28 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 29 ierr = MatSetSizes(A,1,1,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 30 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 31 ierr = MatSetUp(A);CHKERRQ(ierr); 32 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 33 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 34 /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */ 35 ierr = MatShift(A,(PetscReal)1);CHKERRQ(ierr); 36 ierr = MatShift(A,(PetscReal)-1);CHKERRQ(ierr); 37 ierr = TSSetIJacobian(ts,A,A,IJacobian,NULL);CHKERRQ(ierr); 38 ierr = MatDestroy(&A);CHKERRQ(ierr); 39 40 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 41 ierr = VecSetSizes(x,1,PETSC_DECIDE);CHKERRQ(ierr); 42 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 43 ierr = VecSetUp(x);CHKERRQ(ierr); 44 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 45 ierr = VecDestroy(&x);CHKERRQ(ierr); 46 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 47 48 ierr = TSSetStepNumber(ts,0);CHKERRQ(ierr); 49 ierr = TSSetTimeStep(ts,1);CHKERRQ(ierr); 50 ierr = TSSetTime(ts,0);CHKERRQ(ierr); 51 ierr = TSSetMaxTime(ts,PETSC_MAX_REAL);CHKERRQ(ierr); 52 ierr = TSSetMaxSteps(ts,3);CHKERRQ(ierr); 53 54 /* 55 When an ARKIMEX scheme with an explicit stage is used this will error with a message informing the user it is not possible to use 56 a non-trivial mass matrix with ARKIMEX schemes with explicit stages. 57 */ 58 ierr = TSSolve(ts,NULL); 59 if (ierr != PETSC_ERR_ARG_INCOMP) CHKERRQ(ierr); 60 61 ierr = TSDestroy(&ts);CHKERRQ(ierr); 62 ierr = PetscFinalize(); 63 return ierr; 64 } 65 66 PetscErrorCode IFunction(TS ts,PetscReal t,Vec x,Vec xdot,Vec f,void *ctx) 67 { 68 PetscErrorCode ierr; 69 70 PetscFunctionBegin; 71 ierr = VecCopy(xdot,f);CHKERRQ(ierr); 72 ierr = VecScale(f,2.0);CHKERRQ(ierr); 73 ierr = VecShift(f,-1.0);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec x,Vec xdot,PetscReal shift,Mat A,Mat B,void *ctx) 78 { 79 PetscErrorCode ierr; 80 PetscScalar j; 81 82 PetscFunctionBegin; 83 j = shift*2.0; 84 ierr = MatSetValue(B,0,0,j,INSERT_VALUES);CHKERRQ(ierr); 85 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 86 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 /*TEST 91 92 test: 93 suffix: arkimex_explicit_stage 94 requires: defined(PETSC_USE_DEBUG) 95 args: -ts_type arkimex -error_output_stdout 96 filter: egrep -v "(Petsc|on a| at |Configure)" 97 98 test: 99 suffix: arkimex_implicit_stage 100 args: -ts_type arkimex -ts_arkimex_type l2 -ts_monitor_solution -ts_monitor 101 102 TEST*/ 103