Actual source code: ex21.c

      2: static char help[] = "Solves PDE optimization problem.\n\n";

     4:  #include petscda.h
     5:  #include petscpf.h
     6:  #include petscsnes.h

      8: /*

     10:        w - design variables (what we change to get an optimal solution)
     11:        u - state variables (i.e. the PDE solution)
     12:        lambda - the Lagrange multipliers

     14:             U = (w u lambda)

     16:        fu, fw, flambda contain the gradient of L(w,u,lambda)

     18:             FU = (fw fu flambda)

     20:        In this example the PDE is
     21:                              Uxx = 2,
     22:                             u(0) = w(0), thus this is the free parameter
     23:                             u(1) = 0
     24:        the function we wish to minimize is
     25:                             \integral u^{2}

     27:        The exact solution for u is given by u(x) = x*x - 1.25*x + .25

     29:        Use the usual centered finite differences.

     31:        Note we treat the problem as non-linear though it happens to be linear

     33:        See ex22.c for the same code, but that interlaces the u and the lambda

     35: */

     37: typedef struct {
     38:   DA           da1,da2;
     39:   PetscInt     nredundant;
     40:   VecPack      packer;
     41:   PetscViewer  u_viewer,lambda_viewer;
     42:   PetscViewer  fu_viewer,flambda_viewer;
     43: } UserCtx;

     51: int main(int argc,char **argv)
     52: {
     54:   PetscInt       its;
     55:   Vec            U,FU;
     56:   SNES           snes;
     57:   UserCtx        user;

     59:   PetscInitialize(&argc,&argv,(char*)0,help);

     61:   /* Create a global vector that includes a single redundant array and two da arrays */
     62:   VecPackCreate(PETSC_COMM_WORLD,&user.packer);
     63:   user.nredundant = 1;
     64:   VecPackAddArray(user.packer,user.nredundant);
     65:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,1,1,PETSC_NULL,&user.da1);
     66:   VecPackAddDA(user.packer,user.da1);
     67:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-5,1,1,PETSC_NULL,&user.da2);
     68:   VecPackAddDA(user.packer,user.da2);
     69:   VecPackCreateGlobalVector(user.packer,&U);
     70:   VecDuplicate(U,&FU);

     72:   /* create graphics windows */
     73:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);
     74:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"lambda - Lagrange multipliers",-1,-1,-1,-1,&user.lambda_viewer);
     75:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - derivate w.r.t. state variables",-1,-1,-1,-1,&user.fu_viewer);
     76:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"flambda - derivate w.r.t. Lagrange multipliers",-1,-1,-1,-1,&user.flambda_viewer);

     79:   /* create nonlinear solver */
     80:   SNESCreate(PETSC_COMM_WORLD,&snes);
     81:   SNESSetFunction(snes,FU,FormFunction,&user);
     82:   SNESSetFromOptions(snes);
     83:   SNESSetMonitor(snes,Monitor,&user,0);
     84:   SNESSolve(snes,PETSC_NULL,U);
     85:   SNESGetIterationNumber(snes,&its);
     86:   SNESDestroy(snes);

     88:   DADestroy(user.da1);
     89:   DADestroy(user.da2);
     90:   VecPackDestroy(user.packer);
     91:   VecDestroy(U);
     92:   VecDestroy(FU);
     93:   PetscViewerDestroy(user.u_viewer);
     94:   PetscViewerDestroy(user.lambda_viewer);
     95:   PetscViewerDestroy(user.fu_viewer);
     96:   PetscViewerDestroy(user.flambda_viewer);
     97:   PetscFinalize();
     98:   return 0;
     99: }
    100:
    101: /*
    102:       Evaluates FU = Gradiant(L(w,u,lambda))

    104: */
    105: PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void* dummy)
    106: {
    107:   UserCtx        *user = (UserCtx*)dummy;
    109:   PetscInt       xs,xm,i,N;
    110:   PetscScalar    *u,*lambda,*w,*fu,*fw,*flambda,d,h;
    111:   Vec            vu,vlambda,vfu,vflambda;

    114:   VecPackGetLocalVectors(user->packer,&w,&vu,&vlambda);
    115:   VecPackGetLocalVectors(user->packer,&fw,&vfu,&vflambda);
    116:   VecPackScatter(user->packer,U,w,vu,vlambda);

    118:   DAGetCorners(user->da1,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
    119:   DAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0);
    120:   DAVecGetArray(user->da1,vu,&u);
    121:   DAVecGetArray(user->da1,vfu,&fu);
    122:   DAVecGetArray(user->da1,vlambda,&lambda);
    123:   DAVecGetArray(user->da1,vflambda,&flambda);
    124:   d    = (N-1.0);
    125:   h    = 1.0/d;

    127:   /* derivative of L() w.r.t. w */
    128:   if (xs == 0) { /* only first processor computes this */
    129:     fw[0] = -2.*d*lambda[0];
    130:   }

    132:   /* derivative of L() w.r.t. u */
    133:   for (i=xs; i<xs+xm; i++) {
    134:     if      (i == 0)   flambda[0]   =    h*u[0]   + 2.*d*lambda[0]   - d*lambda[1];
    135:     else if (i == 1)   flambda[1]   = 2.*h*u[1]   + 2.*d*lambda[1]   - d*lambda[2];
    136:     else if (i == N-1) flambda[N-1] =    h*u[N-1] + 2.*d*lambda[N-1] - d*lambda[N-2];
    137:     else if (i == N-2) flambda[N-2] = 2.*h*u[N-2] + 2.*d*lambda[N-2] - d*lambda[N-3];
    138:     else               flambda[i]   = 2.*h*u[i]   - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]);
    139:   }

    141:   /* derivative of L() w.r.t. lambda */
    142:   for (i=xs; i<xs+xm; i++) {
    143:     if      (i == 0)   fu[0]   = 2.0*d*(u[0] - w[0]);
    144:     else if (i == N-1) fu[N-1] = 2.0*d*u[N-1];
    145:     else               fu[i]   = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h);
    146:   }

    148:   DAVecRestoreArray(user->da1,vu,&u);
    149:   DAVecRestoreArray(user->da1,vfu,&fu);
    150:   DAVecRestoreArray(user->da1,vlambda,&lambda);
    151:   DAVecRestoreArray(user->da1,vflambda,&flambda);

    153:   VecPackGather(user->packer,FU,fw,vfu,vflambda);
    154:   VecPackRestoreLocalVectors(user->packer,&w,&vu,&vlambda);
    155:   VecPackRestoreLocalVectors(user->packer,&fw,&vfu,&vflambda);
    156:   return(0);
    157: }

    159: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy)
    160: {
    161:   UserCtx        *user = (UserCtx*)dummy;
    163:   PetscScalar    *w;
    164:   Vec            u,lambda,U,F;

    167:   SNESGetSolution(snes,&U);
    168:   VecPackGetAccess(user->packer,U,&w,&u,&lambda);
    169:   VecView(u,user->u_viewer);
    170:   VecView(lambda,user->lambda_viewer);
    171:   VecPackRestoreAccess(user->packer,U,&w,&u,&lambda);

    173:   SNESGetFunction(snes,&F,0,0);
    174:   VecPackGetAccess(user->packer,F,&w,&u,&lambda);
    175:   VecView(u,user->fu_viewer);
    176:   VecView(lambda,user->flambda_viewer);
    177:   VecPackRestoreAccess(user->packer,F,&w,&u,&lambda);
    178:   return(0);
    179: }