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 PetscCall(VecGetArrayRead(X,&x)); 48 PetscCall(VecGetArrayRead(Xdot,&xdot)); 49 PetscCall(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 PetscCall(VecRestoreArrayRead(X,&x)); 53 PetscCall(VecRestoreArrayRead(Xdot,&xdot)); 54 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 70 PetscCall(VecRestoreArrayRead(X,&x)); 71 72 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 73 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 74 if (A != B) { 75 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 76 PetscCall(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 PetscCall(VecGetArrayRead(X,&x)); 91 J[0][0] = 0; 92 J[1][0] = (1.-x[0]*x[0])*x[1]-x[0]; 93 PetscCall(MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES)); 94 PetscCall(VecRestoreArrayRead(X,&x)); 95 96 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 97 PetscCall(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 const PetscScalar *x; 105 PetscReal tfinal, dt; 106 User user = (User)ctx; 107 Vec interpolatedX; 108 109 PetscFunctionBeginUser; 110 PetscCall(TSGetTimeStep(ts,&dt)); 111 PetscCall(TSGetMaxTime(ts,&tfinal)); 112 113 while (user->next_output <= t && user->next_output <= tfinal) { 114 PetscCall(VecDuplicate(X,&interpolatedX)); 115 PetscCall(TSInterpolate(ts,user->next_output,interpolatedX)); 116 PetscCall(VecGetArrayRead(interpolatedX,&x)); 117 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", 118 user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]), 119 (double)PetscRealPart(x[1]))); 120 PetscCall(VecRestoreArrayRead(interpolatedX,&x)); 121 PetscCall(VecDestroy(&interpolatedX)); 122 user->next_output += 0.1; 123 } 124 PetscFunctionReturn(0); 125 } 126 127 int main(int argc,char **argv) 128 { 129 TS ts; 130 PetscBool monitor = PETSC_FALSE; 131 PetscScalar *x_ptr; 132 PetscMPIInt size; 133 struct _n_User user; 134 PetscInt rows,cols; 135 136 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 137 Initialize program 138 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 139 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 140 141 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 142 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 143 144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 145 Set runtime options 146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 147 user.next_output = 0.0; 148 user.mu = 1.0e6; 149 user.steps = 0; 150 user.ftime = 0.5; 151 user.combined = PETSC_FALSE; 152 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 153 PetscCall(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 154 PetscCall(PetscOptionsGetBool(NULL,NULL,"-combined",&user.combined,NULL)); 155 156 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 157 Create necessary matrix and vectors, solve same ODE on every process 158 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 159 rows = 2; 160 cols = user.combined ? 3 : 1; 161 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.Jac)); 162 PetscCall(MatSetSizes(user.Jac,PETSC_DECIDE,PETSC_DECIDE,2,2)); 163 PetscCall(MatSetFromOptions(user.Jac)); 164 PetscCall(MatSetUp(user.Jac)); 165 PetscCall(MatCreateVecs(user.Jac,&user.x,NULL)); 166 167 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 168 Create timestepping solver context 169 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 170 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 171 PetscCall(TSSetType(ts,TSBEULER)); 172 PetscCall(TSSetIFunction(ts,NULL,IFunction,&user)); 173 PetscCall(TSSetIJacobian(ts,user.Jac,user.Jac,IJacobian,&user)); 174 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 175 PetscCall(TSSetMaxTime(ts,user.ftime)); 176 if (monitor) { 177 PetscCall(TSMonitorSet(ts,Monitor,&user,NULL)); 178 } 179 180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 181 Set initial conditions 182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 183 PetscCall(VecGetArray(user.x,&x_ptr)); 184 x_ptr[0] = 2.0; x_ptr[1] = -0.66666654321; 185 PetscCall(VecRestoreArray(user.x,&x_ptr)); 186 PetscCall(TSSetTimeStep(ts,1.0/1024.0)); 187 188 /* Set up forward sensitivity */ 189 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.Jacp)); 190 PetscCall(MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,rows,cols)); 191 PetscCall(MatSetFromOptions(user.Jacp)); 192 PetscCall(MatSetUp(user.Jacp)); 193 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,rows,cols,NULL,&user.sp)); 194 if (user.combined) { 195 PetscCall(MatZeroEntries(user.sp)); 196 PetscCall(MatShift(user.sp,1.0)); 197 } else { 198 PetscCall(MatZeroEntries(user.sp)); 199 } 200 PetscCall(TSForwardSetSensitivities(ts,cols,user.sp)); 201 PetscCall(TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP,&user)); 202 203 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 204 Set runtime options 205 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 206 PetscCall(TSSetFromOptions(ts)); 207 208 PetscCall(TSSolve(ts,user.x)); 209 PetscCall(TSGetSolveTime(ts,&user.ftime)); 210 PetscCall(TSGetStepNumber(ts,&user.steps)); 211 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,user.steps,(double)user.ftime)); 212 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n ode solution \n")); 213 PetscCall(VecView(user.x,PETSC_VIEWER_STDOUT_WORLD)); 214 215 if (user.combined) { 216 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n forward sensitivity: d[y(tf) z(tf)]/d[y0 z0 mu]\n")); 217 } else { 218 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n forward sensitivity: d[y(tf) z(tf)]/d[mu]\n")); 219 } 220 PetscCall(MatView(user.sp,PETSC_VIEWER_STDOUT_WORLD)); 221 222 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 223 Free work space. All PETSc objects should be destroyed when they 224 are no longer needed. 225 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 226 PetscCall(MatDestroy(&user.Jac)); 227 PetscCall(MatDestroy(&user.sp)); 228 PetscCall(MatDestroy(&user.Jacp)); 229 PetscCall(VecDestroy(&user.x)); 230 PetscCall(TSDestroy(&ts)); 231 232 PetscCall(PetscFinalize()); 233 return 0; 234 } 235 236 /*TEST 237 238 test: 239 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -combined 240 requires: !complex !single 241 242 test: 243 suffix: 2 244 requires: !complex !single 245 args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 246 247 TEST*/ 248