# petsc / src / snes / examples / tutorials / ex21.c

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192``` ``` static const char help[] = "Solves PDE optimization problem using full-space method, treats state and adjoint variables separately.\n\n"; #include #include #include #include #include #include /* w - design variables (what we change to get an optimal solution) u - state variables (i.e. the PDE solution) lambda - the Lagrange multipliers U = (w u lambda) fu, fw, flambda contain the gradient of L(w,u,lambda) FU = (fw fu flambda) In this example the PDE is Uxx = 2, u(0) = w(0), thus this is the free parameter u(1) = 0 the function we wish to minimize is \integral u^{2} The exact solution for u is given by u(x) = x*x - 1.25*x + .25 Use the usual centered finite differences. Note we treat the problem as non-linear though it happens to be linear See ex22.c for the same code, but that interlaces the u and the lambda */ typedef struct { DM red1,da1,da2; DM packer; PetscViewer u_viewer,lambda_viewer; PetscViewer fu_viewer,flambda_viewer; } UserCtx; extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*); #undef __FUNCT__ #define __FUNCT__ "main" int main(int argc,char **argv) { PetscErrorCode ierr; PetscInt its; Vec U,FU; SNES snes; UserCtx user; ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr); /* Create a global vector that includes a single redundant array and two da arrays */ ierr = DMCompositeCreate(PETSC_COMM_WORLD,&user.packer);CHKERRQ(ierr); ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1);CHKERRQ(ierr); ierr = DMCompositeAddDM(user.packer,user.red1);CHKERRQ(ierr); ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-5,1,1,NULL,&user.da1);CHKERRQ(ierr); ierr = DMCompositeAddDM(user.packer,user.da1);CHKERRQ(ierr); ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-5,1,1,NULL,&user.da2);CHKERRQ(ierr); ierr = DMCompositeAddDM(user.packer,user.da2);CHKERRQ(ierr); ierr = DMCreateGlobalVector(user.packer,&U);CHKERRQ(ierr); ierr = VecDuplicate(U,&FU);CHKERRQ(ierr); /* create graphics windows */ ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);CHKERRQ(ierr); ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"lambda - Lagrange multipliers",-1,-1,-1,-1,&user.lambda_viewer);CHKERRQ(ierr); ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - derivate w.r.t. state variables",-1,-1,-1,-1,&user.fu_viewer);CHKERRQ(ierr); ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"flambda - derivate w.r.t. Lagrange multipliers",-1,-1,-1,-1,&user.flambda_viewer);CHKERRQ(ierr); /* create nonlinear solver */ ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); ierr = SNESSetFunction(snes,FU,FormFunction,&user);CHKERRQ(ierr); ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); ierr = SNESMonitorSet(snes,Monitor,&user,0);CHKERRQ(ierr); ierr = SNESSolve(snes,NULL,U);CHKERRQ(ierr); ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr); ierr = DMDestroy(&user.red1);CHKERRQ(ierr); ierr = DMDestroy(&user.da1);CHKERRQ(ierr); ierr = DMDestroy(&user.da2);CHKERRQ(ierr); ierr = DMDestroy(&user.packer);CHKERRQ(ierr); ierr = VecDestroy(&U);CHKERRQ(ierr); ierr = VecDestroy(&FU);CHKERRQ(ierr); ierr = PetscViewerDestroy(&user.u_viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&user.lambda_viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&user.fu_viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(&user.flambda_viewer);CHKERRQ(ierr); ierr = PetscFinalize(); return 0; } #undef __FUNCT__ #define __FUNCT__ "FormFunction" /* Evaluates FU = Gradiant(L(w,u,lambda)) */ PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void *dummy) { UserCtx *user = (UserCtx*)dummy; PetscErrorCode ierr; PetscInt xs,xm,i,N; PetscScalar *u,*lambda,*w,*fu,*fw,*flambda,d,h; Vec vw,vu,vlambda,vfw,vfu,vflambda; PetscFunctionBeginUser; ierr = DMCompositeGetLocalVectors(user->packer,&vw,&vu,&vlambda);CHKERRQ(ierr); ierr = DMCompositeGetLocalVectors(user->packer,&vfw,&vfu,&vflambda);CHKERRQ(ierr); ierr = DMCompositeScatter(user->packer,U,vw,vu,vlambda);CHKERRQ(ierr); ierr = DMDAGetCorners(user->da1,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); ierr = DMDAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); ierr = VecGetArray(vw,&w);CHKERRQ(ierr); ierr = VecGetArray(vfw,&fw);CHKERRQ(ierr); ierr = DMDAVecGetArray(user->da1,vu,&u);CHKERRQ(ierr); ierr = DMDAVecGetArray(user->da1,vfu,&fu);CHKERRQ(ierr); ierr = DMDAVecGetArray(user->da1,vlambda,&lambda);CHKERRQ(ierr); ierr = DMDAVecGetArray(user->da1,vflambda,&flambda);CHKERRQ(ierr); d = (N-1.0); h = 1.0/d; /* derivative of L() w.r.t. w */ if (xs == 0) { /* only first processor computes this */ fw[0] = -2.*d*lambda[0]; } /* derivative of L() w.r.t. u */ for (i=xs; ida1,vu,&u);CHKERRQ(ierr); ierr = DMDAVecRestoreArray(user->da1,vfu,&fu);CHKERRQ(ierr); ierr = DMDAVecRestoreArray(user->da1,vlambda,&lambda);CHKERRQ(ierr); ierr = DMDAVecRestoreArray(user->da1,vflambda,&flambda);CHKERRQ(ierr); ierr = DMCompositeGather(user->packer,FU,INSERT_VALUES,vfw,vfu,vflambda);CHKERRQ(ierr); ierr = DMCompositeRestoreLocalVectors(user->packer,&vw,&vu,&vlambda);CHKERRQ(ierr); ierr = DMCompositeRestoreLocalVectors(user->packer,&vfw,&vfu,&vflambda);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "Monitor" PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy) { UserCtx *user = (UserCtx*)dummy; PetscErrorCode ierr; Vec w,u,lambda,U,F; PetscFunctionBeginUser; ierr = SNESGetSolution(snes,&U);CHKERRQ(ierr); ierr = DMCompositeGetAccess(user->packer,U,&w,&u,&lambda);CHKERRQ(ierr); ierr = VecView(u,user->u_viewer);CHKERRQ(ierr); ierr = VecView(lambda,user->lambda_viewer);CHKERRQ(ierr); ierr = DMCompositeRestoreAccess(user->packer,U,&w,&u,&lambda);CHKERRQ(ierr); ierr = SNESGetFunction(snes,&F,0,0);CHKERRQ(ierr); ierr = DMCompositeGetAccess(user->packer,F,&w,&u,&lambda);CHKERRQ(ierr); ierr = VecView(u,user->fu_viewer);CHKERRQ(ierr); ierr = VecView(lambda,user->flambda_viewer);CHKERRQ(ierr); ierr = DMCompositeRestoreAccess(user->packer,F,&w,&u,&lambda);CHKERRQ(ierr); PetscFunctionReturn(0); } ```