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