1 2 static const char help[] = "Solves PDE optimization problem using full-space method, treats state and adjoint variables separately.\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 11 /* 12 13 w - design variables (what we change to get an optimal solution) 14 u - state variables (i.e. the PDE solution) 15 lambda - the Lagrange multipliers 16 17 U = (w u lambda) 18 19 fu, fw, flambda contain the gradient of L(w,u,lambda) 20 21 FU = (fw fu flambda) 22 23 In this example the PDE is 24 Uxx = 2, 25 u(0) = w(0), thus this is the free parameter 26 u(1) = 0 27 the function we wish to minimize is 28 \integral u^{2} 29 30 The exact solution for u is given by u(x) = x*x - 1.25*x + .25 31 32 Use the usual centered finite differences. 33 34 Note we treat the problem as non-linear though it happens to be linear 35 36 See ex22.c for the same code, but that interlaces the u and the lambda 37 38 */ 39 40 typedef struct { 41 DM red1,da1,da2; 42 DM packer; 43 PetscViewer u_viewer,lambda_viewer; 44 PetscViewer fu_viewer,flambda_viewer; 45 } UserCtx; 46 47 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 48 extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*); 49 50 51 int main(int argc,char **argv) 52 { 53 PetscErrorCode ierr; 54 PetscInt its; 55 Vec U,FU; 56 SNES snes; 57 UserCtx user; 58 59 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 60 61 /* Create a global vector that includes a single redundant array and two da arrays */ 62 ierr = DMCompositeCreate(PETSC_COMM_WORLD,&user.packer);CHKERRQ(ierr); 63 ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1);CHKERRQ(ierr); 64 ierr = DMCompositeAddDM(user.packer,user.red1);CHKERRQ(ierr); 65 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,1,1,NULL,&user.da1);CHKERRQ(ierr); 66 ierr = DMSetFromOptions(user.da1);CHKERRQ(ierr); 67 ierr = DMSetUp(user.da1);CHKERRQ(ierr); 68 ierr = DMCompositeAddDM(user.packer,user.da1);CHKERRQ(ierr); 69 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,1,1,NULL,&user.da2);CHKERRQ(ierr); 70 ierr = DMSetFromOptions(user.da2);CHKERRQ(ierr); 71 ierr = DMSetUp(user.da2);CHKERRQ(ierr); 72 ierr = DMDASetFieldName(user.da1,0,"u");CHKERRQ(ierr); 73 ierr = DMDASetFieldName(user.da2,0,"lambda");CHKERRQ(ierr); 74 ierr = DMCompositeAddDM(user.packer,user.da2);CHKERRQ(ierr); 75 ierr = DMCreateGlobalVector(user.packer,&U);CHKERRQ(ierr); 76 ierr = VecDuplicate(U,&FU);CHKERRQ(ierr); 77 78 /* create graphics windows */ 79 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer);CHKERRQ(ierr); 80 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"lambda - Lagrange multipliers",-1,-1,-1,-1,&user.lambda_viewer);CHKERRQ(ierr); 81 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - derivate w.r.t. state variables",-1,-1,-1,-1,&user.fu_viewer);CHKERRQ(ierr); 82 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"flambda - derivate w.r.t. Lagrange multipliers",-1,-1,-1,-1,&user.flambda_viewer);CHKERRQ(ierr); 83 84 85 /* create nonlinear solver */ 86 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 87 ierr = SNESSetFunction(snes,FU,FormFunction,&user);CHKERRQ(ierr); 88 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 89 ierr = SNESMonitorSet(snes,Monitor,&user,0);CHKERRQ(ierr); 90 ierr = SNESSolve(snes,NULL,U);CHKERRQ(ierr); 91 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 92 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 93 94 ierr = DMDestroy(&user.red1);CHKERRQ(ierr); 95 ierr = DMDestroy(&user.da1);CHKERRQ(ierr); 96 ierr = DMDestroy(&user.da2);CHKERRQ(ierr); 97 ierr = DMDestroy(&user.packer);CHKERRQ(ierr); 98 ierr = VecDestroy(&U);CHKERRQ(ierr); 99 ierr = VecDestroy(&FU);CHKERRQ(ierr); 100 ierr = PetscViewerDestroy(&user.u_viewer);CHKERRQ(ierr); 101 ierr = PetscViewerDestroy(&user.lambda_viewer);CHKERRQ(ierr); 102 ierr = PetscViewerDestroy(&user.fu_viewer);CHKERRQ(ierr); 103 ierr = PetscViewerDestroy(&user.flambda_viewer);CHKERRQ(ierr); 104 ierr = PetscFinalize(); 105 return ierr; 106 } 107 108 /* 109 Evaluates FU = Gradiant(L(w,u,lambda)) 110 111 */ 112 PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void *dummy) 113 { 114 UserCtx *user = (UserCtx*)dummy; 115 PetscErrorCode ierr; 116 PetscInt xs,xm,i,N; 117 PetscScalar *u,*lambda,*w,*fu,*fw,*flambda,d,h; 118 Vec vw,vu,vlambda,vfw,vfu,vflambda; 119 120 PetscFunctionBeginUser; 121 ierr = DMCompositeGetLocalVectors(user->packer,&vw,&vu,&vlambda);CHKERRQ(ierr); 122 ierr = DMCompositeGetLocalVectors(user->packer,&vfw,&vfu,&vflambda);CHKERRQ(ierr); 123 ierr = DMCompositeScatter(user->packer,U,vw,vu,vlambda);CHKERRQ(ierr); 124 125 ierr = DMDAGetCorners(user->da1,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 126 ierr = DMDAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 127 ierr = VecGetArray(vw,&w);CHKERRQ(ierr); 128 ierr = VecGetArray(vfw,&fw);CHKERRQ(ierr); 129 ierr = DMDAVecGetArray(user->da1,vu,&u);CHKERRQ(ierr); 130 ierr = DMDAVecGetArray(user->da1,vfu,&fu);CHKERRQ(ierr); 131 ierr = DMDAVecGetArray(user->da1,vlambda,&lambda);CHKERRQ(ierr); 132 ierr = DMDAVecGetArray(user->da1,vflambda,&flambda);CHKERRQ(ierr); 133 d = (N-1.0); 134 h = 1.0/d; 135 136 /* derivative of L() w.r.t. w */ 137 if (xs == 0) { /* only first processor computes this */ 138 fw[0] = -2.*d*lambda[0]; 139 } 140 141 /* derivative of L() w.r.t. u */ 142 for (i=xs; i<xs+xm; i++) { 143 if (i == 0) flambda[0] = h*u[0] + 2.*d*lambda[0] - d*lambda[1]; 144 else if (i == 1) flambda[1] = 2.*h*u[1] + 2.*d*lambda[1] - d*lambda[2]; 145 else if (i == N-1) flambda[N-1] = h*u[N-1] + 2.*d*lambda[N-1] - d*lambda[N-2]; 146 else if (i == N-2) flambda[N-2] = 2.*h*u[N-2] + 2.*d*lambda[N-2] - d*lambda[N-3]; 147 else flambda[i] = 2.*h*u[i] - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]); 148 } 149 150 /* derivative of L() w.r.t. lambda */ 151 for (i=xs; i<xs+xm; i++) { 152 if (i == 0) fu[0] = 2.0*d*(u[0] - w[0]); 153 else if (i == N-1) fu[N-1] = 2.0*d*u[N-1]; 154 else fu[i] = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h); 155 } 156 157 ierr = VecRestoreArray(vw,&w);CHKERRQ(ierr); 158 ierr = VecRestoreArray(vfw,&fw);CHKERRQ(ierr); 159 ierr = DMDAVecRestoreArray(user->da1,vu,&u);CHKERRQ(ierr); 160 ierr = DMDAVecRestoreArray(user->da1,vfu,&fu);CHKERRQ(ierr); 161 ierr = DMDAVecRestoreArray(user->da1,vlambda,&lambda);CHKERRQ(ierr); 162 ierr = DMDAVecRestoreArray(user->da1,vflambda,&flambda);CHKERRQ(ierr); 163 164 ierr = DMCompositeGather(user->packer,INSERT_VALUES,FU,vfw,vfu,vflambda);CHKERRQ(ierr); 165 ierr = DMCompositeRestoreLocalVectors(user->packer,&vw,&vu,&vlambda);CHKERRQ(ierr); 166 ierr = DMCompositeRestoreLocalVectors(user->packer,&vfw,&vfu,&vflambda);CHKERRQ(ierr); 167 PetscFunctionReturn(0); 168 } 169 170 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy) 171 { 172 UserCtx *user = (UserCtx*)dummy; 173 PetscErrorCode ierr; 174 Vec w,u,lambda,U,F; 175 176 PetscFunctionBeginUser; 177 ierr = SNESGetSolution(snes,&U);CHKERRQ(ierr); 178 ierr = DMCompositeGetAccess(user->packer,U,&w,&u,&lambda);CHKERRQ(ierr); 179 ierr = VecView(u,user->u_viewer);CHKERRQ(ierr); 180 ierr = VecView(lambda,user->lambda_viewer);CHKERRQ(ierr); 181 ierr = DMCompositeRestoreAccess(user->packer,U,&w,&u,&lambda);CHKERRQ(ierr); 182 183 ierr = SNESGetFunction(snes,&F,0,0);CHKERRQ(ierr); 184 ierr = DMCompositeGetAccess(user->packer,F,&w,&u,&lambda);CHKERRQ(ierr); 185 ierr = VecView(u,user->fu_viewer);CHKERRQ(ierr); 186 ierr = VecView(lambda,user->flambda_viewer);CHKERRQ(ierr); 187 ierr = DMCompositeRestoreAccess(user->packer,F,&w,&u,&lambda);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 /*TEST 192 193 test: 194 nsize: 4 195 args: -snes_linesearch_monitor -snes_mf -pc_type none -snes_monitor_short -nox -ksp_monitor_short -snes_converged_reason 196 requires: !single 197 198 TEST*/ 199