Actual source code: ex26.c

  1: static char help[] = "Solves the trival ODE 2 du/dt = 1, u(0) = 0. \n\n";

  3: #include <petscts.h>
  4: #include <petscpc.h>

  6: PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*);
  7: PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

  9: int main(int argc,char **argv)
 10: {
 11:   TS              ts;
 12:   Vec             x;
 13:   Vec             f;
 14:   Mat             A;
 15:   PetscErrorCode  ierr;

 17:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

 19:   TSCreate(PETSC_COMM_WORLD,&ts);
 20:   TSSetEquationType(ts,TS_EQ_ODE_IMPLICIT);
 21:   VecCreate(PETSC_COMM_WORLD,&f);
 22:   VecSetSizes(f,1,PETSC_DECIDE);
 23:   VecSetFromOptions(f);
 24:   VecSetUp(f);
 25:   TSSetIFunction(ts,f,IFunction,NULL);
 26:   VecDestroy(&f);

 28:   MatCreate(PETSC_COMM_WORLD,&A);
 29:   MatSetSizes(A,1,1,PETSC_DECIDE,PETSC_DECIDE);
 30:   MatSetFromOptions(A);
 31:   MatSetUp(A);
 32:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 33:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 34:   /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
 35:   MatShift(A,(PetscReal)1);
 36:   MatShift(A,(PetscReal)-1);
 37:   TSSetIJacobian(ts,A,A,IJacobian,NULL);
 38:   MatDestroy(&A);

 40:   VecCreate(PETSC_COMM_WORLD,&x);
 41:   VecSetSizes(x,1,PETSC_DECIDE);
 42:   VecSetFromOptions(x);
 43:   VecSetUp(x);
 44:   TSSetSolution(ts,x);
 45:   VecDestroy(&x);
 46:   TSSetFromOptions(ts);

 48:   TSSetStepNumber(ts,0);
 49:   TSSetTimeStep(ts,1);
 50:   TSSetTime(ts,0);
 51:   TSSetMaxTime(ts,PETSC_MAX_REAL);
 52:   TSSetMaxSteps(ts,3);

 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:   TSSolve(ts,NULL);
 59:   if (ierr != PETSC_ERR_ARG_INCOMP) 

 61:   TSDestroy(&ts);
 62:   PetscFinalize();
 63:   return ierr;
 64: }

 66: PetscErrorCode IFunction(TS ts,PetscReal t,Vec x,Vec xdot,Vec f,void *ctx)
 67: {

 71:   VecCopy(xdot,f);
 72:   VecScale(f,2.0);
 73:   VecShift(f,-1.0);
 74:   return(0);
 75: }

 77: PetscErrorCode IJacobian(TS ts,PetscReal t,Vec x,Vec xdot,PetscReal shift,Mat A,Mat B,void *ctx)
 78: {
 80:   PetscScalar    j;

 83:   j = shift*2.0;
 84:   MatSetValue(B,0,0,j,INSERT_VALUES);
 85:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 86:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 87:   return(0);
 88: }


 91: /*TEST

 93:     test:
 94:       suffix: arkimex_explicit_stage
 95:       requires: define(PETSC_USE_DEBUG)
 96:       args: -ts_type arkimex -error_output_stdout
 97:       filter:  egrep -v "(Petsc|on a| at |Configure)"

 99:     test:
100:       suffix: arkimex_implicit_stage
101:       args: -ts_type arkimex -ts_arkimex_type l2 -ts_monitor_solution -ts_monitor

103: TEST*/