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 UserCtx user; 77 DM red,da; 78 SNES snes; 79 DM packer; 80 PetscBool use_monitor = PETSC_FALSE; 81 82 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 83 84 /* Hardwire several options; can be changed at command line */ 85 PetscCall(PetscOptionsInsertString(NULL,common_options)); 86 PetscCall(PetscOptionsInsertString(NULL,matrix_free_options)); 87 PetscCall(PetscOptionsInsert(NULL,&argc,&argv,NULL)); 88 PetscCall(PetscOptionsGetBool(NULL,NULL,"-use_monitor",&use_monitor,PETSC_IGNORE)); 89 90 /* Create a global vector that includes a single redundant array and two da arrays */ 91 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD,&packer)); 92 PetscCall(DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red)); 93 PetscCall(DMSetOptionsPrefix(red,"red_")); 94 PetscCall(DMCompositeAddDM(packer,red)); 95 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,2,1,NULL,&da)); 96 PetscCall(DMSetOptionsPrefix(red,"da_")); 97 PetscCall(DMSetFromOptions(da)); 98 PetscCall(DMSetUp(da)); 99 PetscCall(DMDASetFieldName(da,0,"u")); 100 PetscCall(DMDASetFieldName(da,1,"lambda")); 101 PetscCall(DMCompositeAddDM(packer,(DM)da)); 102 PetscCall(DMSetApplicationContext(packer,&user)); 103 104 packer->ops->creatematrix = DMCreateMatrix_MF; 105 106 /* create nonlinear multi-level solver */ 107 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 108 PetscCall(SNESSetDM(snes,packer)); 109 PetscCall(SNESSetFunction(snes,NULL,ComputeFunction,NULL)); 110 PetscCall(SNESSetJacobian(snes,NULL, NULL,ComputeJacobian_MF,NULL)); 111 112 PetscCall(SNESSetFromOptions(snes)); 113 114 if (use_monitor) { 115 /* create graphics windows */ 116 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer)); 117 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu_lambda - derivate w.r.t. state variables and Lagrange multipliers",-1,-1,-1,-1,&user.fu_lambda_viewer)); 118 PetscCall(SNESMonitorSet(snes,Monitor,0,0)); 119 } 120 121 PetscCall(SNESSolve(snes,NULL,NULL)); 122 PetscCall(SNESDestroy(&snes)); 123 124 PetscCall(DMDestroy(&red)); 125 PetscCall(DMDestroy(&da)); 126 PetscCall(DMDestroy(&packer)); 127 if (use_monitor) { 128 PetscCall(PetscViewerDestroy(&user.u_lambda_viewer)); 129 PetscCall(PetscViewerDestroy(&user.fu_lambda_viewer)); 130 } 131 PetscCall(PetscFinalize()); 132 return 0; 133 } 134 135 typedef struct { 136 PetscScalar u; 137 PetscScalar lambda; 138 } ULambda; 139 140 /* 141 Evaluates FU = Gradiant(L(w,u,lambda)) 142 143 This local function acts on the ghosted version of U (accessed via DMCompositeGetLocalVectors() and 144 DMCompositeScatter()) BUT the global, nonghosted version of FU (via DMCompositeGetAccess()). 145 146 */ 147 PetscErrorCode ComputeFunction(SNES snes,Vec U,Vec FU,void *ctx) 148 { 149 PetscInt xs,xm,i,N; 150 ULambda *u_lambda,*fu_lambda; 151 PetscScalar d,h,*w,*fw; 152 Vec vw,vfw,vu_lambda,vfu_lambda; 153 DM packer,red,da; 154 155 PetscFunctionBeginUser; 156 PetscCall(VecGetDM(U, &packer)); 157 PetscCall(DMCompositeGetEntries(packer,&red,&da)); 158 PetscCall(DMCompositeGetLocalVectors(packer,&vw,&vu_lambda)); 159 PetscCall(DMCompositeScatter(packer,U,vw,vu_lambda)); 160 PetscCall(DMCompositeGetAccess(packer,FU,&vfw,&vfu_lambda)); 161 162 PetscCall(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 163 PetscCall(DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0)); 164 PetscCall(VecGetArray(vw,&w)); 165 PetscCall(VecGetArray(vfw,&fw)); 166 PetscCall(DMDAVecGetArray(da,vu_lambda,&u_lambda)); 167 PetscCall(DMDAVecGetArray(da,vfu_lambda,&fu_lambda)); 168 d = N-1.0; 169 h = 1.0/d; 170 171 /* derivative of L() w.r.t. w */ 172 if (xs == 0) { /* only first processor computes this */ 173 fw[0] = -2.0*d*u_lambda[0].lambda; 174 } 175 176 /* derivative of L() w.r.t. u */ 177 for (i=xs; i<xs+xm; i++) { 178 if (i == 0) fu_lambda[0].lambda = h*u_lambda[0].u + 2.*d*u_lambda[0].lambda - d*u_lambda[1].lambda; 179 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; 180 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; 181 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; 182 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); 183 } 184 185 /* derivative of L() w.r.t. lambda */ 186 for (i=xs; i<xs+xm; i++) { 187 if (i == 0) fu_lambda[0].u = 2.0*d*(u_lambda[0].u - w[0]); 188 else if (i == N-1) fu_lambda[N-1].u = 2.0*d*u_lambda[N-1].u; 189 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); 190 } 191 192 PetscCall(VecRestoreArray(vw,&w)); 193 PetscCall(VecRestoreArray(vfw,&fw)); 194 PetscCall(DMDAVecRestoreArray(da,vu_lambda,&u_lambda)); 195 PetscCall(DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda)); 196 PetscCall(DMCompositeRestoreLocalVectors(packer,&vw,&vu_lambda)); 197 PetscCall(DMCompositeRestoreAccess(packer,FU,&vfw,&vfu_lambda)); 198 PetscCall(PetscLogFlops(13.0*N)); 199 PetscFunctionReturn(0); 200 } 201 202 /* 203 Computes the exact solution 204 */ 205 PetscErrorCode u_solution(void *dummy,PetscInt n,const PetscScalar *x,PetscScalar *u) 206 { 207 PetscInt i; 208 209 PetscFunctionBeginUser; 210 for (i=0; i<n; i++) u[2*i] = x[i]*x[i] - 1.25*x[i] + .25; 211 PetscFunctionReturn(0); 212 } 213 214 PetscErrorCode ExactSolution(DM packer,Vec U) 215 { 216 PF pf; 217 Vec x,u_global; 218 PetscScalar *w; 219 DM da; 220 PetscInt m; 221 222 PetscFunctionBeginUser; 223 PetscCall(DMCompositeGetEntries(packer,&m,&da)); 224 225 PetscCall(PFCreate(PETSC_COMM_WORLD,1,2,&pf)); 226 /* The cast through PETSC_UINTPTR_T is so that compilers will warn about casting to void * from void(*)(void) */ 227 PetscCall(PFSetType(pf,PFQUICK,(void*)(PETSC_UINTPTR_T)u_solution)); 228 PetscCall(DMGetCoordinates(da,&x)); 229 if (!x) { 230 PetscCall(DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0)); 231 PetscCall(DMGetCoordinates(da,&x)); 232 } 233 PetscCall(DMCompositeGetAccess(packer,U,&w,&u_global,0)); 234 if (w) w[0] = .25; 235 PetscCall(PFApplyVec(pf,x,u_global)); 236 PetscCall(PFDestroy(&pf)); 237 PetscCall(DMCompositeRestoreAccess(packer,U,&w,&u_global,0)); 238 PetscFunctionReturn(0); 239 } 240 241 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy) 242 { 243 UserCtx *user; 244 PetscInt m,N; 245 PetscScalar *w,*dw; 246 Vec u_lambda,U,F,Uexact; 247 DM packer; 248 PetscReal norm; 249 DM da; 250 251 PetscFunctionBeginUser; 252 PetscCall(SNESGetDM(snes,&packer)); 253 PetscCall(DMGetApplicationContext(packer,&user)); 254 PetscCall(SNESGetSolution(snes,&U)); 255 PetscCall(DMCompositeGetAccess(packer,U,&w,&u_lambda)); 256 PetscCall(VecView(u_lambda,user->u_lambda_viewer)); 257 PetscCall(DMCompositeRestoreAccess(packer,U,&w,&u_lambda)); 258 259 PetscCall(SNESGetFunction(snes,&F,0,0)); 260 PetscCall(DMCompositeGetAccess(packer,F,&w,&u_lambda)); 261 /* ierr = VecView(u_lambda,user->fu_lambda_viewer); */ 262 PetscCall(DMCompositeRestoreAccess(packer,U,&w,&u_lambda)); 263 264 PetscCall(DMCompositeGetEntries(packer,&m,&da)); 265 PetscCall(DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0)); 266 PetscCall(VecDuplicate(U,&Uexact)); 267 PetscCall(ExactSolution(packer,Uexact)); 268 PetscCall(VecAXPY(Uexact,-1.0,U)); 269 PetscCall(DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda)); 270 PetscCall(VecStrideNorm(u_lambda,0,NORM_2,&norm)); 271 norm = norm/PetscSqrtReal((PetscReal)N-1.); 272 if (dw) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Error at x = 0 %g\n",(double)norm,(double)PetscRealPart(dw[0]))); 273 PetscCall(VecView(u_lambda,user->fu_lambda_viewer)); 274 PetscCall(DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda)); 275 PetscCall(VecDestroy(&Uexact)); 276 PetscFunctionReturn(0); 277 } 278 279 PetscErrorCode DMCreateMatrix_MF(DM packer,Mat *A) 280 { 281 Vec t; 282 PetscInt m; 283 284 PetscFunctionBeginUser; 285 PetscCall(DMGetGlobalVector(packer,&t)); 286 PetscCall(VecGetLocalSize(t,&m)); 287 PetscCall(DMRestoreGlobalVector(packer,&t)); 288 PetscCall(MatCreateMFFD(PETSC_COMM_WORLD,m,m,PETSC_DETERMINE,PETSC_DETERMINE,A)); 289 PetscCall(MatSetUp(*A)); 290 PetscCall(MatSetDM(*A,packer)); 291 PetscFunctionReturn(0); 292 } 293 294 PetscErrorCode ComputeJacobian_MF(SNES snes,Vec x,Mat A,Mat B,void *ctx) 295 { 296 PetscFunctionBeginUser; 297 PetscCall(MatMFFDSetFunction(A,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes)); 298 PetscCall(MatMFFDSetBase(A,x,NULL)); 299 PetscFunctionReturn(0); 300 } 301 302 /*TEST 303 304 test: 305 nsize: 2 306 args: -da_grid_x 10 -snes_converged_reason -ksp_converged_reason -snes_view 307 requires: !single 308 309 TEST*/ 310