Actual source code: ex3sa.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      2: static char help[] = "Adjoint and tangent linear sensitivity analysis of the basic equation for generator stability analysis.\n";

\begin{eqnarray}
\frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
\frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
\end{eqnarray}

     13: /*
     14:   This code demonstrate the sensitivity analysis interface to a system of ordinary differential equations with discontinuities.
     15:   It computes the sensitivities of an integral cost function
     16:   \int c*max(0,\theta(t)-u_s)^beta dt
     17:   w.r.t. initial conditions and the parameter P_m.
     18:   Backward Euler method is used for time integration.
     19:   The discontinuities are detected with TSEvent.
     20:  */

     22: #include <petscts.h>
     23: #include "ex3.h"

     25: int main(int argc,char **argv)
     26: {
     27:   TS             ts,quadts;     /* ODE integrator */
     28:   Vec            U;             /* solution will be stored here */
     30:   PetscMPIInt    size;
     31:   PetscInt       n = 2;
     32:   AppCtx         ctx;
     33:   PetscScalar    *u;
     34:   PetscReal      du[2] = {0.0,0.0};
     35:   PetscBool      ensemble = PETSC_FALSE,flg1,flg2;
     36:   PetscReal      ftime;
     37:   PetscInt       steps;
     38:   PetscScalar    *x_ptr,*y_ptr,*s_ptr;
     39:   Vec            lambda[1],q,mu[1];
     40:   PetscInt       direction[2];
     41:   PetscBool      terminate[2];
     42:   Mat            qgrad;
     43:   Mat            sp;            /* Forward sensitivity matrix */
     44:   SAMethod       sa;

     46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     47:      Initialize program
     48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     49:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     50:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     51:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");

     53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     54:     Create necessary matrix and vectors
     55:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     56:   MatCreate(PETSC_COMM_WORLD,&ctx.Jac);
     57:   MatSetSizes(ctx.Jac,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
     58:   MatSetType(ctx.Jac,MATDENSE);
     59:   MatSetFromOptions(ctx.Jac);
     60:   MatSetUp(ctx.Jac);
     61:   MatCreateVecs(ctx.Jac,&U,NULL);
     62:   MatCreate(PETSC_COMM_WORLD,&ctx.Jacp);
     63:   MatSetSizes(ctx.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
     64:   MatSetFromOptions(ctx.Jacp);
     65:   MatSetUp(ctx.Jacp);
     66:   MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&ctx.DRDP);
     67:   MatSetUp(ctx.DRDP);
     68:   MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&ctx.DRDU);
     69:   MatSetUp(ctx.DRDU);

     71:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     72:     Set runtime options
     73:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     74:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
     75:   {
     76:     ctx.beta    = 2;
     77:     ctx.c       = 10000.0;
     78:     ctx.u_s     = 1.0;
     79:     ctx.omega_s = 1.0;
     80:     ctx.omega_b = 120.0*PETSC_PI;
     81:     ctx.H       = 5.0;
     82:     PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
     83:     ctx.D       = 5.0;
     84:     PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
     85:     ctx.E       = 1.1378;
     86:     ctx.V       = 1.0;
     87:     ctx.X       = 0.545;
     88:     ctx.Pmax    = ctx.E*ctx.V/ctx.X;
     89:     ctx.Pmax_ini = ctx.Pmax;
     90:     PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
     91:     ctx.Pm      = 1.1;
     92:     PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
     93:     ctx.tf      = 0.1;
     94:     ctx.tcl     = 0.2;
     95:     PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
     96:     PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
     97:     PetscOptionsBool("-ensemble","Run ensemble of different initial conditions","",ensemble,&ensemble,NULL);
     98:     if (ensemble) {
     99:       ctx.tf      = -1;
    100:       ctx.tcl     = -1;
    101:     }

    103:     VecGetArray(U,&u);
    104:     u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
    105:     u[1] = 1.0;
    106:     PetscOptionsRealArray("-u","Initial solution","",u,&n,&flg1);
    107:     n    = 2;
    108:     PetscOptionsRealArray("-du","Perturbation in initial solution","",du,&n,&flg2);
    109:     u[0] += du[0];
    110:     u[1] += du[1];
    111:     VecRestoreArray(U,&u);
    112:     if (flg1 || flg2) {
    113:       ctx.tf      = -1;
    114:       ctx.tcl     = -1;
    115:     }
    116:     sa = SA_ADJ;
    117:     PetscOptionsEnum("-sa_method","Sensitivity analysis method (adj or tlm)","",SAMethods,(PetscEnum)sa,(PetscEnum*)&sa,NULL);
    118:   }
    119:   PetscOptionsEnd();

    121:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    122:      Create timestepping solver context
    123:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    124:   TSCreate(PETSC_COMM_WORLD,&ts);
    125:   TSSetProblemType(ts,TS_NONLINEAR);
    126:   TSSetType(ts,TSBEULER);
    127:   TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
    128:   TSSetRHSJacobian(ts,ctx.Jac,ctx.Jac,(TSRHSJacobian)RHSJacobian,&ctx);

    130:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    131:      Set initial conditions
    132:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    133:   TSSetSolution(ts,U);

    135:   /*   Set RHS JacobianP */
    136:   TSSetRHSJacobianP(ts,ctx.Jacp,RHSJacobianP,&ctx);

    138:   TSCreateQuadratureTS(ts,PETSC_FALSE,&quadts);
    139:   TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);
    140:   TSSetRHSJacobian(quadts,ctx.DRDU,ctx.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);
    141:   TSSetRHSJacobianP(quadts,ctx.DRDP,DRDPJacobianTranspose,&ctx);
    142:   if (sa == SA_ADJ) {
    143:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    144:       Save trajectory of solution so that TSAdjointSolve() may be used
    145:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    146:     TSSetSaveTrajectory(ts);
    147:     MatCreateVecs(ctx.Jac,&lambda[0],NULL);
    148:     MatCreateVecs(ctx.Jacp,&mu[0],NULL);
    149:     TSSetCostGradients(ts,1,lambda,mu);
    150:   }

    152:   if (sa == SA_TLM) {
    153:     PetscScalar val[2];
    154:     PetscInt    row[]={0,1},col[]={0};

    156:     MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&qgrad);
    157:     MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&sp);
    158:     TSForwardSetSensitivities(ts,1,sp);
    159:     TSForwardSetSensitivities(quadts,1,qgrad);
    160:     val[0] = 1./PetscSqrtScalar(1.-(ctx.Pm/ctx.Pmax)*(ctx.Pm/ctx.Pmax))/ctx.Pmax;
    161:     val[1] = 0.0;
    162:     MatSetValues(sp,2,row,1,col,val,INSERT_VALUES);
    163:     MatAssemblyBegin(sp,MAT_FINAL_ASSEMBLY);
    164:     MatAssemblyEnd(sp,MAT_FINAL_ASSEMBLY);
    165:   }

    167:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    168:      Set solver options
    169:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    170:   TSSetMaxTime(ts,1.0);
    171:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
    172:   TSSetTimeStep(ts,0.03125);
    173:   TSSetFromOptions(ts);

    175:   direction[0] = direction[1] = 1;
    176:   terminate[0] = terminate[1] = PETSC_FALSE;

    178:   TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&ctx);

    180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    181:      Solve nonlinear system
    182:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    183:   if (ensemble) {
    184:     for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
    185:       VecGetArray(U,&u);
    186:       u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
    187:       u[1] = ctx.omega_s;
    188:       u[0] += du[0];
    189:       u[1] += du[1];
    190:       VecRestoreArray(U,&u);
    191:       TSSetTimeStep(ts,0.03125);
    192:       TSSolve(ts,U);
    193:     }
    194:   } else {
    195:     TSSolve(ts,U);
    196:   }
    197:   TSGetSolveTime(ts,&ftime);
    198:   TSGetStepNumber(ts,&steps);

    200:   if (sa == SA_ADJ) {
    201:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    202:        Adjoint model starts here
    203:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    204:     /*   Set initial conditions for the adjoint integration */
    205:     VecGetArray(lambda[0],&y_ptr);
    206:     y_ptr[0] = 0.0; y_ptr[1] = 0.0;
    207:     VecRestoreArray(lambda[0],&y_ptr);

    209:     VecGetArray(mu[0],&x_ptr);
    210:     x_ptr[0] = 0.0;
    211:     VecRestoreArray(mu[0],&x_ptr);

    213:     TSAdjointSolve(ts);

    215:     PetscPrintf(PETSC_COMM_WORLD,"\n lambda: d[Psi(tf)]/d[phi0]  d[Psi(tf)]/d[omega0]\n");
    216:     VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD);
    217:     PetscPrintf(PETSC_COMM_WORLD,"\n mu: d[Psi(tf)]/d[pm]\n");
    218:     VecView(mu[0],PETSC_VIEWER_STDOUT_WORLD);
    219:     TSGetCostIntegral(ts,&q);
    220:     VecGetArray(q,&x_ptr);
    221:     PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(x_ptr[0]-ctx.Pm));
    222:     VecRestoreArray(q,&x_ptr);
    223:     ComputeSensiP(lambda[0],mu[0],&ctx);
    224:     VecGetArray(mu[0],&x_ptr);
    225:     PetscPrintf(PETSC_COMM_WORLD,"\n gradient=%g\n",(double)x_ptr[0]);
    226:     VecRestoreArray(mu[0],&x_ptr);
    227:     VecDestroy(&lambda[0]);
    228:     VecDestroy(&mu[0]);
    229:   }
    230:   if (sa == SA_TLM) {
    231:     PetscPrintf(PETSC_COMM_WORLD,"\n trajectory sensitivity: d[phi(tf)]/d[pm]  d[omega(tf)]/d[pm]\n");
    232:     MatView(sp,PETSC_VIEWER_STDOUT_WORLD);
    233:     TSGetCostIntegral(ts,&q);
    234:     VecGetArray(q,&s_ptr);
    235:     PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(s_ptr[0]-ctx.Pm));
    236:     VecRestoreArray(q,&s_ptr);
    237:     MatDenseGetArray(qgrad,&s_ptr);
    238:     PetscPrintf(PETSC_COMM_WORLD,"\n gradient=%g\n",(double)s_ptr[0]);
    239:     MatDenseRestoreArray(qgrad,&s_ptr);
    240:     MatDestroy(&qgrad);
    241:     MatDestroy(&sp);
    242:   }
    243:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    244:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
    245:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    246:   MatDestroy(&ctx.Jac);
    247:   MatDestroy(&ctx.Jacp);
    248:   MatDestroy(&ctx.DRDU);
    249:   MatDestroy(&ctx.DRDP);
    250:   VecDestroy(&U);
    251:   TSDestroy(&ts);
    252:   PetscFinalize();
    253:   return ierr;
    254: }

    257: /*TEST

    259:    build:
    260:       requires: !complex !single

    262:    test:
    263:       args: -sa_method adj -viewer_binary_skip_info -ts_type cn -pc_type lu

    265:    test:
    266:       suffix: 2
    267:       args: -sa_method tlm -ts_type cn -pc_type lu

    269:    test:
    270:       suffix: 3
    271:       args: -sa_method adj -ts_type rk -ts_rk_type 2a -ts_adapt_type dsp

    273: TEST*/