1 2 static const char help[] = "Solves PDE optimization problem using full-space method, interlaces state and adjoint variables.\n\n"; 3 4 5 #include <petscdm.h> 6 #include <petscdmda.h> 7 #include <petscdmredundant.h> 8 #include <petscdmcomposite.h> 9 #include <petscpf.h> 10 #include <petscsnes.h> 11 #include <petsc/private/dmimpl.h> 12 13 /* 14 15 w - design variables (what we change to get an optimal solution) 16 u - state variables (i.e. the PDE solution) 17 lambda - the Lagrange multipliers 18 19 U = (w [u_0 lambda_0 u_1 lambda_1 .....]) 20 21 fu, fw, flambda contain the gradient of L(w,u,lambda) 22 23 FU = (fw [fu_0 flambda_0 .....]) 24 25 In this example the PDE is 26 Uxx = 2, 27 u(0) = w(0), thus this is the free parameter 28 u(1) = 0 29 the function we wish to minimize is 30 \integral u^{2} 31 32 The exact solution for u is given by u(x) = x*x - 1.25*x + .25 33 34 Use the usual centered finite differences. 35 36 Note we treat the problem as non-linear though it happens to be linear 37 38 See ex21.c for the same code, but that does NOT interlaces the u and the lambda 39 40 The vectors u_lambda and fu_lambda contain the u and the lambda interlaced 41 */ 42 43 typedef struct { 44 PetscViewer u_lambda_viewer; 45 PetscViewer fu_lambda_viewer; 46 } UserCtx; 47 48 extern PetscErrorCode ComputeFunction(SNES,Vec,Vec,void*); 49 extern PetscErrorCode ComputeJacobian_MF(SNES,Vec,Mat,Mat,void*); 50 extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*); 51 52 /* 53 Uses full multigrid preconditioner with GMRES (with no preconditioner inside the GMRES) as the 54 smoother on all levels. This is because (1) in the matrix free case no matrix entries are 55 available for doing Jacobi or SOR preconditioning and (2) the explicit matrix case the diagonal 56 entry for the control variable is zero which means default SOR will not work. 57 58 */ 59 char common_options[] = "-ksp_type fgmres\ 60 -snes_grid_sequence 2 \ 61 -pc_type mg\ 62 -mg_levels_pc_type none \ 63 -mg_coarse_pc_type none \ 64 -pc_mg_type full \ 65 -mg_coarse_ksp_type gmres \ 66 -mg_levels_ksp_type gmres \ 67 -mg_coarse_ksp_max_it 6 \ 68 -mg_levels_ksp_max_it 3"; 69 70 char matrix_free_options[] = "-mat_mffd_compute_normu no \ 71 -mat_mffd_type wp"; 72 73 extern PetscErrorCode DMCreateMatrix_MF(DM,Mat*); 74 75 int main(int argc,char **argv) 76 { 77 PetscErrorCode ierr; 78 UserCtx user; 79 DM red,da; 80 SNES snes; 81 DM packer; 82 PetscBool use_monitor = PETSC_FALSE; 83 84 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 85 ierr = PetscOptionsSetFromOptions(NULL);CHKERRQ(ierr); 86 87 /* Hardwire several options; can be changed at command line */ 88 ierr = PetscOptionsInsertString(NULL,common_options);CHKERRQ(ierr); 89 ierr = PetscOptionsInsertString(NULL,matrix_free_options);CHKERRQ(ierr); 90 ierr = PetscOptionsInsert(NULL,&argc,&argv,NULL);CHKERRQ(ierr); 91 ierr = PetscOptionsGetBool(NULL,NULL,"-use_monitor",&use_monitor,PETSC_IGNORE);CHKERRQ(ierr); 92 93 /* Create a global vector that includes a single redundant array and two da arrays */ 94 ierr = DMCompositeCreate(PETSC_COMM_WORLD,&packer);CHKERRQ(ierr); 95 ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red);CHKERRQ(ierr); 96 ierr = DMSetOptionsPrefix(red,"red_");CHKERRQ(ierr); 97 ierr = DMCompositeAddDM(packer,red);CHKERRQ(ierr); 98 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,2,1,NULL,&da);CHKERRQ(ierr); 99 ierr = DMSetOptionsPrefix(red,"da_");CHKERRQ(ierr); 100 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 101 ierr = DMSetUp(da);CHKERRQ(ierr); 102 ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr); 103 ierr = DMDASetFieldName(da,1,"lambda");CHKERRQ(ierr); 104 ierr = DMCompositeAddDM(packer,(DM)da);CHKERRQ(ierr); 105 ierr = DMSetApplicationContext(packer,&user);CHKERRQ(ierr); 106 107 packer->ops->creatematrix = DMCreateMatrix_MF; 108 109 /* create nonlinear multi-level solver */ 110 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 111 ierr = SNESSetDM(snes,packer);CHKERRQ(ierr); 112 ierr = SNESSetFunction(snes,NULL,ComputeFunction,NULL);CHKERRQ(ierr); 113 ierr = SNESSetJacobian(snes,NULL, NULL,ComputeJacobian_MF,NULL);CHKERRQ(ierr); 114 115 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 116 117 if (use_monitor) { 118 /* create graphics windows */ 119 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer);CHKERRQ(ierr); 120 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer);CHKERRQ(ierr); 121 ierr = SNESMonitorSet(snes,Monitor,0,0);CHKERRQ(ierr); 122 } 123 124 ierr = SNESSolve(snes,NULL,NULL);CHKERRQ(ierr); 125 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 126 127 ierr = DMDestroy(&red);CHKERRQ(ierr); 128 ierr = DMDestroy(&da);CHKERRQ(ierr); 129 ierr = DMDestroy(&packer);CHKERRQ(ierr); 130 if (use_monitor) { 131 ierr = PetscViewerDestroy(&user.u_lambda_viewer);CHKERRQ(ierr); 132 ierr = PetscViewerDestroy(&user.fu_lambda_viewer);CHKERRQ(ierr); 133 } 134 ierr = PetscFinalize(); 135 return ierr; 136 } 137 138 typedef struct { 139 PetscScalar u; 140 PetscScalar lambda; 141 } ULambda; 142 143 /* 144 Evaluates FU = Gradiant(L(w,u,lambda)) 145 146 This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and 147 DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()). 148 149 */ 150 PetscErrorCode ComputeFunction(SNES snes,Vec U,Vec FU,void *ctx) 151 { 152 PetscErrorCode ierr; 153 PetscInt xs,xm,i,N; 154 ULambda *u_lambda,*fu_lambda; 155 PetscScalar d,h,*w,*fw; 156 Vec vw,vfw,vu_lambda,vfu_lambda; 157 DM packer,red,da; 158 159 PetscFunctionBeginUser; 160 ierr = VecGetDM(U, &packer);CHKERRQ(ierr); 161 ierr = DMCompositeGetEntries(packer,&red,&da);CHKERRQ(ierr); 162 ierr = DMCompositeGetLocalVectors(packer,&vw,&vu_lambda);CHKERRQ(ierr); 163 ierr = DMCompositeScatter(packer,U,vw,vu_lambda);CHKERRQ(ierr); 164 ierr = DMCompositeGetAccess(packer,FU,&vfw,&vfu_lambda);CHKERRQ(ierr); 165 166 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 167 ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 168 ierr = VecGetArray(vw,&w);CHKERRQ(ierr); 169 ierr = VecGetArray(vfw,&fw);CHKERRQ(ierr); 170 ierr = DMDAVecGetArray(da,vu_lambda,&u_lambda);CHKERRQ(ierr); 171 ierr = DMDAVecGetArray(da,vfu_lambda,&fu_lambda);CHKERRQ(ierr); 172 d = N-1.0; 173 h = 1.0/d; 174 175 /* derivative of L() w.r.t. w */ 176 if (xs == 0) { /* only first processor computes this */ 177 fw[0] = -2.0*d*u_lambda[0].lambda; 178 } 179 180 /* derivative of L() w.r.t. u */ 181 for (i=xs; i<xs+xm; i++) { 182 if (i == 0) fu_lambda[0].lambda = h*u_lambda[0].u + 2.*d*u_lambda[0].lambda - d*u_lambda[1].lambda; 183 else if (i == 1) fu_lambda[1].lambda = 2.*h*u_lambda[1].u + 2.*d*u_lambda[1].lambda - d*u_lambda[2].lambda; 184 else if (i == N-1) fu_lambda[N-1].lambda = h*u_lambda[N-1].u + 2.*d*u_lambda[N-1].lambda - d*u_lambda[N-2].lambda; 185 else if (i == N-2) fu_lambda[N-2].lambda = 2.*h*u_lambda[N-2].u + 2.*d*u_lambda[N-2].lambda - d*u_lambda[N-3].lambda; 186 else fu_lambda[i].lambda = 2.*h*u_lambda[i].u - d*(u_lambda[i+1].lambda - 2.0*u_lambda[i].lambda + u_lambda[i-1].lambda); 187 } 188 189 /* derivative of L() w.r.t. lambda */ 190 for (i=xs; i<xs+xm; i++) { 191 if (i == 0) fu_lambda[0].u = 2.0*d*(u_lambda[0].u - w[0]); 192 else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u; 193 else fu_lambda[i].u = -(d*(u_lambda[i+1].u - 2.0*u_lambda[i].u + u_lambda[i-1].u) - 2.0*h); 194 } 195 196 ierr = VecRestoreArray(vw,&w);CHKERRQ(ierr); 197 ierr = VecRestoreArray(vfw,&fw);CHKERRQ(ierr); 198 ierr = DMDAVecRestoreArray(da,vu_lambda,&u_lambda);CHKERRQ(ierr); 199 ierr = DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda);CHKERRQ(ierr); 200 ierr = DMCompositeRestoreLocalVectors(packer,&vw,&vu_lambda);CHKERRQ(ierr); 201 ierr = DMCompositeRestoreAccess(packer,FU,&vfw,&vfu_lambda);CHKERRQ(ierr); 202 ierr = PetscLogFlops(13.0*N);CHKERRQ(ierr); 203 PetscFunctionReturn(0); 204 } 205 206 /* 207 Computes the exact solution 208 */ 209 PetscErrorCode u_solution(void *dummy,PetscInt n,const PetscScalar *x,PetscScalar *u) 210 { 211 PetscInt i; 212 213 PetscFunctionBeginUser; 214 for (i=0; i<n; i++) u[2*i] = x[i]*x[i] - 1.25*x[i] + .25; 215 PetscFunctionReturn(0); 216 } 217 218 PetscErrorCode ExactSolution(DM packer,Vec U) 219 { 220 PF pf; 221 Vec x,u_global; 222 PetscScalar *w; 223 DM da; 224 PetscErrorCode ierr; 225 PetscInt m; 226 227 PetscFunctionBeginUser; 228 ierr = DMCompositeGetEntries(packer,&m,&da);CHKERRQ(ierr); 229 230 ierr = PFCreate(PETSC_COMM_WORLD,1,2,&pf);CHKERRQ(ierr); 231 /* The cast through PETSC_UINTPTR_T is so that compilers will warn about casting to void * from void(*)(void) */ 232 ierr = PFSetType(pf,PFQUICK,(void*)(PETSC_UINTPTR_T)u_solution);CHKERRQ(ierr); 233 ierr = DMGetCoordinates(da,&x);CHKERRQ(ierr); 234 if (!x) { 235 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr); 236 ierr = DMGetCoordinates(da,&x);CHKERRQ(ierr); 237 } 238 ierr = DMCompositeGetAccess(packer,U,&w,&u_global,0);CHKERRQ(ierr); 239 if (w) w[0] = .25; 240 ierr = PFApplyVec(pf,x,u_global);CHKERRQ(ierr); 241 ierr = PFDestroy(&pf);CHKERRQ(ierr); 242 ierr = DMCompositeRestoreAccess(packer,U,&w,&u_global,0);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy) 247 { 248 UserCtx *user; 249 PetscErrorCode ierr; 250 PetscInt m,N; 251 PetscScalar *w,*dw; 252 Vec u_lambda,U,F,Uexact; 253 DM packer; 254 PetscReal norm; 255 DM da; 256 257 PetscFunctionBeginUser; 258 ierr = SNESGetDM(snes,&packer);CHKERRQ(ierr); 259 ierr = DMGetApplicationContext(packer,&user);CHKERRQ(ierr); 260 ierr = SNESGetSolution(snes,&U);CHKERRQ(ierr); 261 ierr = DMCompositeGetAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr); 262 ierr = VecView(u_lambda,user->u_lambda_viewer);CHKERRQ(ierr); 263 ierr = DMCompositeRestoreAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr); 264 265 ierr = SNESGetFunction(snes,&F,0,0);CHKERRQ(ierr); 266 ierr = DMCompositeGetAccess(packer,F,&w,&u_lambda);CHKERRQ(ierr); 267 /* ierr = VecView(u_lambda,user->fu_lambda_viewer); */ 268 ierr = DMCompositeRestoreAccess(packer,U,&w,&u_lambda);CHKERRQ(ierr); 269 270 ierr = DMCompositeGetEntries(packer,&m,&da);CHKERRQ(ierr); 271 ierr = DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 272 ierr = VecDuplicate(U,&Uexact);CHKERRQ(ierr); 273 ierr = ExactSolution(packer,Uexact);CHKERRQ(ierr); 274 ierr = VecAXPY(Uexact,-1.0,U);CHKERRQ(ierr); 275 ierr = DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda);CHKERRQ(ierr); 276 ierr = VecStrideNorm(u_lambda,0,NORM_2,&norm);CHKERRQ(ierr); 277 norm = norm/PetscSqrtReal((PetscReal)N-1.); 278 if (dw) ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Error at x = 0 %g\n",(double)norm,(double)PetscRealPart(dw[0]));CHKERRQ(ierr); 279 ierr = VecView(u_lambda,user->fu_lambda_viewer);CHKERRQ(ierr); 280 ierr = DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda);CHKERRQ(ierr); 281 ierr = VecDestroy(&Uexact);CHKERRQ(ierr); 282 PetscFunctionReturn(0); 283 } 284 285 PetscErrorCode DMCreateMatrix_MF(DM packer,Mat *A) 286 { 287 PetscErrorCode ierr; 288 Vec t; 289 PetscInt m; 290 291 PetscFunctionBeginUser; 292 ierr = DMGetGlobalVector(packer,&t);CHKERRQ(ierr); 293 ierr = VecGetLocalSize(t,&m);CHKERRQ(ierr); 294 ierr = DMRestoreGlobalVector(packer,&t);CHKERRQ(ierr); 295 ierr = MatCreateMFFD(PETSC_COMM_WORLD,m,m,PETSC_DETERMINE,PETSC_DETERMINE,A);CHKERRQ(ierr); 296 ierr = MatSetUp(*A);CHKERRQ(ierr); 297 ierr = MatSetDM(*A,packer);CHKERRQ(ierr); 298 PetscFunctionReturn(0); 299 } 300 301 PetscErrorCode ComputeJacobian_MF(SNES snes,Vec x,Mat A,Mat B,void *ctx) 302 { 303 PetscErrorCode ierr; 304 305 PetscFunctionBeginUser; 306 ierr = MatMFFDSetFunction(A,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes);CHKERRQ(ierr); 307 ierr = MatMFFDSetBase(A,x,NULL);CHKERRQ(ierr); 308 PetscFunctionReturn(0); 309 } 310 311 312 /*TEST 313 314 test: 315 nsize: 2 316 args: -da_grid_x 10 -snes_converged_reason -ksp_converged_reason -snes_view 317 requires: !single 318 319 TEST*/ 320