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