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