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 const PetscScalar *x; 105 PetscReal tfinal, dt; 106 User user = (User)ctx; 107 Vec interpolatedX; 108 109 PetscFunctionBeginUser; 110 CHKERRQ(TSGetTimeStep(ts,&dt)); 111 CHKERRQ(TSGetMaxTime(ts,&tfinal)); 112 113 while (user->next_output <= t && user->next_output <= tfinal) { 114 CHKERRQ(VecDuplicate(X,&interpolatedX)); 115 CHKERRQ(TSInterpolate(ts,user->next_output,interpolatedX)); 116 CHKERRQ(VecGetArrayRead(interpolatedX,&x)); 117 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(interpolatedX,&x)); 121 CHKERRQ(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 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 140 141 CHKERRMPI(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 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 153 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 154 CHKERRQ(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 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.Jac)); 162 CHKERRQ(MatSetSizes(user.Jac,PETSC_DECIDE,PETSC_DECIDE,2,2)); 163 CHKERRQ(MatSetFromOptions(user.Jac)); 164 CHKERRQ(MatSetUp(user.Jac)); 165 CHKERRQ(MatCreateVecs(user.Jac,&user.x,NULL)); 166 167 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 168 Create timestepping solver context 169 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 170 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 171 CHKERRQ(TSSetType(ts,TSBEULER)); 172 CHKERRQ(TSSetIFunction(ts,NULL,IFunction,&user)); 173 CHKERRQ(TSSetIJacobian(ts,user.Jac,user.Jac,IJacobian,&user)); 174 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 175 CHKERRQ(TSSetMaxTime(ts,user.ftime)); 176 if (monitor) { 177 CHKERRQ(TSMonitorSet(ts,Monitor,&user,NULL)); 178 } 179 180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 181 Set initial conditions 182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 183 CHKERRQ(VecGetArray(user.x,&x_ptr)); 184 x_ptr[0] = 2.0; x_ptr[1] = -0.66666654321; 185 CHKERRQ(VecRestoreArray(user.x,&x_ptr)); 186 CHKERRQ(TSSetTimeStep(ts,1.0/1024.0)); 187 188 /* Set up forward sensitivity */ 189 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.Jacp)); 190 CHKERRQ(MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,rows,cols)); 191 CHKERRQ(MatSetFromOptions(user.Jacp)); 192 CHKERRQ(MatSetUp(user.Jacp)); 193 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,rows,cols,NULL,&user.sp)); 194 if (user.combined) { 195 CHKERRQ(MatZeroEntries(user.sp)); 196 CHKERRQ(MatShift(user.sp,1.0)); 197 } else { 198 CHKERRQ(MatZeroEntries(user.sp)); 199 } 200 CHKERRQ(TSForwardSetSensitivities(ts,cols,user.sp)); 201 CHKERRQ(TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP,&user)); 202 203 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 204 Set runtime options 205 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 206 CHKERRQ(TSSetFromOptions(ts)); 207 208 CHKERRQ(TSSolve(ts,user.x)); 209 CHKERRQ(TSGetSolveTime(ts,&user.ftime)); 210 CHKERRQ(TSGetStepNumber(ts,&user.steps)); 211 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,user.steps,(double)user.ftime)); 212 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n ode solution \n")); 213 CHKERRQ(VecView(user.x,PETSC_VIEWER_STDOUT_WORLD)); 214 215 if (user.combined) { 216 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n forward sensitivity: d[y(tf) z(tf)]/d[y0 z0 mu]\n")); 217 } else { 218 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n forward sensitivity: d[y(tf) z(tf)]/d[mu]\n")); 219 } 220 CHKERRQ(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 CHKERRQ(MatDestroy(&user.Jac)); 227 CHKERRQ(MatDestroy(&user.sp)); 228 CHKERRQ(MatDestroy(&user.Jacp)); 229 CHKERRQ(VecDestroy(&user.x)); 230 CHKERRQ(TSDestroy(&ts)); 231 232 CHKERRQ(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