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 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 83 84 /* Hardwire several options; can be changed at command line */ 85 CHKERRQ(PetscOptionsInsertString(NULL,common_options)); 86 CHKERRQ(PetscOptionsInsertString(NULL,matrix_free_options)); 87 CHKERRQ(PetscOptionsInsert(NULL,&argc,&argv,NULL)); 88 CHKERRQ(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 CHKERRQ(DMCompositeCreate(PETSC_COMM_WORLD,&packer)); 92 CHKERRQ(DMRedundantCreate(PETSC_COMM_WORLD,0,1,&red)); 93 CHKERRQ(DMSetOptionsPrefix(red,"red_")); 94 CHKERRQ(DMCompositeAddDM(packer,red)); 95 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,2,1,NULL,&da)); 96 CHKERRQ(DMSetOptionsPrefix(red,"da_")); 97 CHKERRQ(DMSetFromOptions(da)); 98 CHKERRQ(DMSetUp(da)); 99 CHKERRQ(DMDASetFieldName(da,0,"u")); 100 CHKERRQ(DMDASetFieldName(da,1,"lambda")); 101 CHKERRQ(DMCompositeAddDM(packer,(DM)da)); 102 CHKERRQ(DMSetApplicationContext(packer,&user)); 103 104 packer->ops->creatematrix = DMCreateMatrix_MF; 105 106 /* create nonlinear multi-level solver */ 107 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 108 CHKERRQ(SNESSetDM(snes,packer)); 109 CHKERRQ(SNESSetFunction(snes,NULL,ComputeFunction,NULL)); 110 CHKERRQ(SNESSetJacobian(snes,NULL, NULL,ComputeJacobian_MF,NULL)); 111 112 CHKERRQ(SNESSetFromOptions(snes)); 113 114 if (use_monitor) { 115 /* create graphics windows */ 116 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u_lambda - state variables and Lagrange multipliers",-1,-1,-1,-1,&user.u_lambda_viewer)); 117 CHKERRQ(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 CHKERRQ(SNESMonitorSet(snes,Monitor,0,0)); 119 } 120 121 CHKERRQ(SNESSolve(snes,NULL,NULL)); 122 CHKERRQ(SNESDestroy(&snes)); 123 124 CHKERRQ(DMDestroy(&red)); 125 CHKERRQ(DMDestroy(&da)); 126 CHKERRQ(DMDestroy(&packer)); 127 if (use_monitor) { 128 CHKERRQ(PetscViewerDestroy(&user.u_lambda_viewer)); 129 CHKERRQ(PetscViewerDestroy(&user.fu_lambda_viewer)); 130 } 131 CHKERRQ(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 CHKERRQ(VecGetDM(U, &packer)); 157 CHKERRQ(DMCompositeGetEntries(packer,&red,&da)); 158 CHKERRQ(DMCompositeGetLocalVectors(packer,&vw,&vu_lambda)); 159 CHKERRQ(DMCompositeScatter(packer,U,vw,vu_lambda)); 160 CHKERRQ(DMCompositeGetAccess(packer,FU,&vfw,&vfu_lambda)); 161 162 CHKERRQ(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 163 CHKERRQ(DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0)); 164 CHKERRQ(VecGetArray(vw,&w)); 165 CHKERRQ(VecGetArray(vfw,&fw)); 166 CHKERRQ(DMDAVecGetArray(da,vu_lambda,&u_lambda)); 167 CHKERRQ(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 CHKERRQ(VecRestoreArray(vw,&w)); 193 CHKERRQ(VecRestoreArray(vfw,&fw)); 194 CHKERRQ(DMDAVecRestoreArray(da,vu_lambda,&u_lambda)); 195 CHKERRQ(DMDAVecRestoreArray(da,vfu_lambda,&fu_lambda)); 196 CHKERRQ(DMCompositeRestoreLocalVectors(packer,&vw,&vu_lambda)); 197 CHKERRQ(DMCompositeRestoreAccess(packer,FU,&vfw,&vfu_lambda)); 198 CHKERRQ(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 CHKERRQ(DMCompositeGetEntries(packer,&m,&da)); 224 225 CHKERRQ(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 CHKERRQ(PFSetType(pf,PFQUICK,(void*)(PETSC_UINTPTR_T)u_solution)); 228 CHKERRQ(DMGetCoordinates(da,&x)); 229 if (!x) { 230 CHKERRQ(DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0)); 231 CHKERRQ(DMGetCoordinates(da,&x)); 232 } 233 CHKERRQ(DMCompositeGetAccess(packer,U,&w,&u_global,0)); 234 if (w) w[0] = .25; 235 CHKERRQ(PFApplyVec(pf,x,u_global)); 236 CHKERRQ(PFDestroy(&pf)); 237 CHKERRQ(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 CHKERRQ(SNESGetDM(snes,&packer)); 253 CHKERRQ(DMGetApplicationContext(packer,&user)); 254 CHKERRQ(SNESGetSolution(snes,&U)); 255 CHKERRQ(DMCompositeGetAccess(packer,U,&w,&u_lambda)); 256 CHKERRQ(VecView(u_lambda,user->u_lambda_viewer)); 257 CHKERRQ(DMCompositeRestoreAccess(packer,U,&w,&u_lambda)); 258 259 CHKERRQ(SNESGetFunction(snes,&F,0,0)); 260 CHKERRQ(DMCompositeGetAccess(packer,F,&w,&u_lambda)); 261 /* ierr = VecView(u_lambda,user->fu_lambda_viewer); */ 262 CHKERRQ(DMCompositeRestoreAccess(packer,U,&w,&u_lambda)); 263 264 CHKERRQ(DMCompositeGetEntries(packer,&m,&da)); 265 CHKERRQ(DMDAGetInfo(da,0,&N,0,0,0,0,0,0,0,0,0,0,0)); 266 CHKERRQ(VecDuplicate(U,&Uexact)); 267 CHKERRQ(ExactSolution(packer,Uexact)); 268 CHKERRQ(VecAXPY(Uexact,-1.0,U)); 269 CHKERRQ(DMCompositeGetAccess(packer,Uexact,&dw,&u_lambda)); 270 CHKERRQ(VecStrideNorm(u_lambda,0,NORM_2,&norm)); 271 norm = norm/PetscSqrtReal((PetscReal)N-1.); 272 if (dw) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Error at x = 0 %g\n",(double)norm,(double)PetscRealPart(dw[0]))); 273 CHKERRQ(VecView(u_lambda,user->fu_lambda_viewer)); 274 CHKERRQ(DMCompositeRestoreAccess(packer,Uexact,&dw,&u_lambda)); 275 CHKERRQ(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 CHKERRQ(DMGetGlobalVector(packer,&t)); 286 CHKERRQ(VecGetLocalSize(t,&m)); 287 CHKERRQ(DMRestoreGlobalVector(packer,&t)); 288 CHKERRQ(MatCreateMFFD(PETSC_COMM_WORLD,m,m,PETSC_DETERMINE,PETSC_DETERMINE,A)); 289 CHKERRQ(MatSetUp(*A)); 290 CHKERRQ(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 CHKERRQ(MatMFFDSetFunction(A,(PetscErrorCode (*)(void*,Vec,Vec))SNESComputeFunction,snes)); 298 CHKERRQ(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