1 #define c11 1.0 2 #define c12 0 3 #define c21 2.0 4 #define c22 1.0 5 static char help[] = "Solves the van der Pol equation.\n\ 6 Input parameters include:\n"; 7 8 /* 9 Concepts: TS^forward sensitivity analysis for time-dependent nonlinear problems 10 Concepts: TS^van der Pol equation DAE equivalent 11 Processors: 1 12 */ 13 /* ------------------------------------------------------------------------ 14 15 This code demonstrates how to compute trajectory sensitivties w.r.t. the stiffness parameter mu. 16 1) Use two vectors s and sp for sensitivities w.r.t. initial values and paraeters respectively. This case requires the original Jacobian matrix and a JacobianP matrix for the only parameter mu. 17 2) Consider the initial values to be parameters as well. Then there are three parameters in total. The JacobianP matrix will be combined matrix of the Jacobian matrix and JacobianP matrix in the previous case. This choice can be selected by using command line option '-combined' 18 19 ------------------------------------------------------------------------- */ 20 #include <petscts.h> 21 #include <petsctao.h> 22 23 typedef struct _n_User *User; 24 struct _n_User { 25 PetscReal mu; 26 PetscReal next_output; 27 PetscBool combined; 28 /* Sensitivity analysis support */ 29 PetscInt steps; 30 PetscReal ftime; 31 Mat Jac; /* Jacobian matrix */ 32 Mat Jacp; /* JacobianP matrix */ 33 Vec x; 34 Mat sp; /* forward sensitivity variables */ 35 }; 36 37 /* 38 User-defined routines 39 */ 40 static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 41 { 42 User user = (User)ctx; 43 const PetscScalar *x,*xdot; 44 PetscScalar *f; 45 46 PetscFunctionBeginUser; 47 CHKERRQ(VecGetArrayRead(X,&x)); 48 CHKERRQ(VecGetArrayRead(Xdot,&xdot)); 49 CHKERRQ(VecGetArray(F,&f)); 50 f[0] = xdot[0] - x[1]; 51 f[1] = c21*(xdot[0]-x[1]) + xdot[1] - user->mu*((1.0-x[0]*x[0])*x[1] - x[0]) ; 52 CHKERRQ(VecRestoreArrayRead(X,&x)); 53 CHKERRQ(VecRestoreArrayRead(Xdot,&xdot)); 54 CHKERRQ(VecRestoreArray(F,&f)); 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) 59 { 60 User user = (User)ctx; 61 PetscInt rowcol[] = {0,1}; 62 PetscScalar J[2][2]; 63 const PetscScalar *x; 64 65 PetscFunctionBeginUser; 66 CHKERRQ(VecGetArrayRead(X,&x)); 67 J[0][0] = a; J[0][1] = -1.0; 68 J[1][0] = c21*a + user->mu*(1.0 + 2.0*x[0]*x[1]); J[1][1] = -c21 + a - user->mu*(1.0-x[0]*x[0]); 69 CHKERRQ(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 70 CHKERRQ(VecRestoreArrayRead(X,&x)); 71 72 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 73 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 74 if (A != B) { 75 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 76 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 77 } 78 PetscFunctionReturn(0); 79 } 80 81 static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx) 82 { 83 User user = (User)ctx; 84 PetscInt row[] = {0,1},col[]={0}; 85 PetscScalar J[2][1]; 86 const PetscScalar *x; 87 88 PetscFunctionBeginUser; 89 if (user->combined) col[0] = 2; 90 CHKERRQ(VecGetArrayRead(X,&x)); 91 J[0][0] = 0; 92 J[1][0] = (1.-x[0]*x[0])*x[1]-x[0]; 93 CHKERRQ(MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES)); 94 CHKERRQ(VecRestoreArrayRead(X,&x)); 95 96 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 97 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 98 PetscFunctionReturn(0); 99 } 100 101 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 102 static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx) 103 { 104 PetscErrorCode ierr; 105 const PetscScalar *x; 106 PetscReal tfinal, dt; 107 User user = (User)ctx; 108 Vec interpolatedX; 109 110 PetscFunctionBeginUser; 111 CHKERRQ(TSGetTimeStep(ts,&dt)); 112 CHKERRQ(TSGetMaxTime(ts,&tfinal)); 113 114 while (user->next_output <= t && user->next_output <= tfinal) { 115 CHKERRQ(VecDuplicate(X,&interpolatedX)); 116 CHKERRQ(TSInterpolate(ts,user->next_output,interpolatedX)); 117 CHKERRQ(VecGetArrayRead(interpolatedX,&x)); 118 ierr = PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", 119 user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]), 120 (double)PetscRealPart(x[1]));CHKERRQ(ierr); 121 CHKERRQ(VecRestoreArrayRead(interpolatedX,&x)); 122 CHKERRQ(VecDestroy(&interpolatedX)); 123 user->next_output += 0.1; 124 } 125 PetscFunctionReturn(0); 126 } 127 128 int main(int argc,char **argv) 129 { 130 TS ts; 131 PetscBool monitor = PETSC_FALSE; 132 PetscScalar *x_ptr; 133 PetscMPIInt size; 134 struct _n_User user; 135 PetscInt rows,cols; 136 PetscErrorCode ierr; 137 138 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139 Initialize program 140 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 141 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 142 143 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 144 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 145 146 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147 Set runtime options 148 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 149 user.next_output = 0.0; 150 user.mu = 1.0e6; 151 user.steps = 0; 152 user.ftime = 0.5; 153 user.combined = PETSC_FALSE; 154 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 155 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 156 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-combined",&user.combined,NULL)); 157 158 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 159 Create necessary matrix and vectors, solve same ODE on every process 160 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 161 rows = 2; 162 cols = user.combined ? 3 : 1; 163 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.Jac)); 164 CHKERRQ(MatSetSizes(user.Jac,PETSC_DECIDE,PETSC_DECIDE,2,2)); 165 CHKERRQ(MatSetFromOptions(user.Jac)); 166 CHKERRQ(MatSetUp(user.Jac)); 167 CHKERRQ(MatCreateVecs(user.Jac,&user.x,NULL)); 168 169 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 170 Create timestepping solver context 171 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 172 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 173 CHKERRQ(TSSetType(ts,TSBEULER)); 174 CHKERRQ(TSSetIFunction(ts,NULL,IFunction,&user)); 175 CHKERRQ(TSSetIJacobian(ts,user.Jac,user.Jac,IJacobian,&user)); 176 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 177 CHKERRQ(TSSetMaxTime(ts,user.ftime)); 178 if (monitor) { 179 CHKERRQ(TSMonitorSet(ts,Monitor,&user,NULL)); 180 } 181 182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 183 Set initial conditions 184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 185 CHKERRQ(VecGetArray(user.x,&x_ptr)); 186 x_ptr[0] = 2.0; x_ptr[1] = -0.66666654321; 187 CHKERRQ(VecRestoreArray(user.x,&x_ptr)); 188 CHKERRQ(TSSetTimeStep(ts,1.0/1024.0)); 189 190 /* Set up forward sensitivity */ 191 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.Jacp)); 192 CHKERRQ(MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,rows,cols)); 193 CHKERRQ(MatSetFromOptions(user.Jacp)); 194 CHKERRQ(MatSetUp(user.Jacp)); 195 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,rows,cols,NULL,&user.sp)); 196 if (user.combined) { 197 CHKERRQ(MatZeroEntries(user.sp)); 198 CHKERRQ(MatShift(user.sp,1.0)); 199 } else { 200 CHKERRQ(MatZeroEntries(user.sp)); 201 } 202 CHKERRQ(TSForwardSetSensitivities(ts,cols,user.sp)); 203 CHKERRQ(TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP,&user)); 204 205 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 206 Set runtime options 207 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 208 CHKERRQ(TSSetFromOptions(ts)); 209 210 CHKERRQ(TSSolve(ts,user.x)); 211 CHKERRQ(TSGetSolveTime(ts,&user.ftime)); 212 CHKERRQ(TSGetStepNumber(ts,&user.steps)); 213 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,user.steps,(double)user.ftime)); 214 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n ode solution \n")); 215 CHKERRQ(VecView(user.x,PETSC_VIEWER_STDOUT_WORLD)); 216 217 if (user.combined) { 218 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n forward sensitivity: d[y(tf) z(tf)]/d[y0 z0 mu]\n")); 219 } else { 220 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n forward sensitivity: d[y(tf) z(tf)]/d[mu]\n")); 221 } 222 CHKERRQ(MatView(user.sp,PETSC_VIEWER_STDOUT_WORLD)); 223 224 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 225 Free work space. All PETSc objects should be destroyed when they 226 are no longer needed. 227 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 228 CHKERRQ(MatDestroy(&user.Jac)); 229 CHKERRQ(MatDestroy(&user.sp)); 230 CHKERRQ(MatDestroy(&user.Jacp)); 231 CHKERRQ(VecDestroy(&user.x)); 232 CHKERRQ(TSDestroy(&ts)); 233 234 ierr = PetscFinalize(); 235 return(ierr); 236 } 237 238 /*TEST 239 240 test: 241 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -combined 242 requires: !complex !single 243 244 test: 245 suffix: 2 246 requires: !complex !single 247 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 248 249 TEST*/ 250