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