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 int main(int argc,char **argv) 51 { 52 PetscInt its; 53 Vec U,FU; 54 SNES snes; 55 UserCtx user; 56 57 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 58 59 /* Create a global vector that includes a single redundant array and two da arrays */ 60 CHKERRQ(DMCompositeCreate(PETSC_COMM_WORLD,&user.packer)); 61 CHKERRQ(DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1)); 62 CHKERRQ(DMCompositeAddDM(user.packer,user.red1)); 63 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,1,1,NULL,&user.da1)); 64 CHKERRQ(DMSetFromOptions(user.da1)); 65 CHKERRQ(DMSetUp(user.da1)); 66 CHKERRQ(DMCompositeAddDM(user.packer,user.da1)); 67 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,5,1,1,NULL,&user.da2)); 68 CHKERRQ(DMSetFromOptions(user.da2)); 69 CHKERRQ(DMSetUp(user.da2)); 70 CHKERRQ(DMDASetFieldName(user.da1,0,"u")); 71 CHKERRQ(DMDASetFieldName(user.da2,0,"lambda")); 72 CHKERRQ(DMCompositeAddDM(user.packer,user.da2)); 73 CHKERRQ(DMCreateGlobalVector(user.packer,&U)); 74 CHKERRQ(VecDuplicate(U,&FU)); 75 76 /* create graphics windows */ 77 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"u - state variables",-1,-1,-1,-1,&user.u_viewer)); 78 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"lambda - Lagrange multipliers",-1,-1,-1,-1,&user.lambda_viewer)); 79 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"fu - derivate w.r.t. state variables",-1,-1,-1,-1,&user.fu_viewer)); 80 CHKERRQ(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"flambda - derivate w.r.t. Lagrange multipliers",-1,-1,-1,-1,&user.flambda_viewer)); 81 82 /* create nonlinear solver */ 83 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 84 CHKERRQ(SNESSetFunction(snes,FU,FormFunction,&user)); 85 CHKERRQ(SNESSetFromOptions(snes)); 86 CHKERRQ(SNESMonitorSet(snes,Monitor,&user,0)); 87 CHKERRQ(SNESSolve(snes,NULL,U)); 88 CHKERRQ(SNESGetIterationNumber(snes,&its)); 89 CHKERRQ(SNESDestroy(&snes)); 90 91 CHKERRQ(DMDestroy(&user.red1)); 92 CHKERRQ(DMDestroy(&user.da1)); 93 CHKERRQ(DMDestroy(&user.da2)); 94 CHKERRQ(DMDestroy(&user.packer)); 95 CHKERRQ(VecDestroy(&U)); 96 CHKERRQ(VecDestroy(&FU)); 97 CHKERRQ(PetscViewerDestroy(&user.u_viewer)); 98 CHKERRQ(PetscViewerDestroy(&user.lambda_viewer)); 99 CHKERRQ(PetscViewerDestroy(&user.fu_viewer)); 100 CHKERRQ(PetscViewerDestroy(&user.flambda_viewer)); 101 CHKERRQ(PetscFinalize()); 102 return 0; 103 } 104 105 /* 106 Evaluates FU = Gradiant(L(w,u,lambda)) 107 108 */ 109 PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void *dummy) 110 { 111 UserCtx *user = (UserCtx*)dummy; 112 PetscInt xs,xm,i,N; 113 PetscScalar *u,*lambda,*w,*fu,*fw,*flambda,d,h; 114 Vec vw,vu,vlambda,vfw,vfu,vflambda; 115 116 PetscFunctionBeginUser; 117 CHKERRQ(DMCompositeGetLocalVectors(user->packer,&vw,&vu,&vlambda)); 118 CHKERRQ(DMCompositeGetLocalVectors(user->packer,&vfw,&vfu,&vflambda)); 119 CHKERRQ(DMCompositeScatter(user->packer,U,vw,vu,vlambda)); 120 121 CHKERRQ(DMDAGetCorners(user->da1,&xs,NULL,NULL,&xm,NULL,NULL)); 122 CHKERRQ(DMDAGetInfo(user->da1,0,&N,0,0,0,0,0,0,0,0,0,0,0)); 123 CHKERRQ(VecGetArray(vw,&w)); 124 CHKERRQ(VecGetArray(vfw,&fw)); 125 CHKERRQ(DMDAVecGetArray(user->da1,vu,&u)); 126 CHKERRQ(DMDAVecGetArray(user->da1,vfu,&fu)); 127 CHKERRQ(DMDAVecGetArray(user->da1,vlambda,&lambda)); 128 CHKERRQ(DMDAVecGetArray(user->da1,vflambda,&flambda)); 129 d = (N-1.0); 130 h = 1.0/d; 131 132 /* derivative of L() w.r.t. w */ 133 if (xs == 0) { /* only first processor computes this */ 134 fw[0] = -2.*d*lambda[0]; 135 } 136 137 /* derivative of L() w.r.t. u */ 138 for (i=xs; i<xs+xm; i++) { 139 if (i == 0) flambda[0] = h*u[0] + 2.*d*lambda[0] - d*lambda[1]; 140 else if (i == 1) flambda[1] = 2.*h*u[1] + 2.*d*lambda[1] - d*lambda[2]; 141 else if (i == N-1) flambda[N-1] = h*u[N-1] + 2.*d*lambda[N-1] - d*lambda[N-2]; 142 else if (i == N-2) flambda[N-2] = 2.*h*u[N-2] + 2.*d*lambda[N-2] - d*lambda[N-3]; 143 else flambda[i] = 2.*h*u[i] - d*(lambda[i+1] - 2.0*lambda[i] + lambda[i-1]); 144 } 145 146 /* derivative of L() w.r.t. lambda */ 147 for (i=xs; i<xs+xm; i++) { 148 if (i == 0) fu[0] = 2.0*d*(u[0] - w[0]); 149 else if (i == N-1) fu[N-1] = 2.0*d*u[N-1]; 150 else fu[i] = -(d*(u[i+1] - 2.0*u[i] + u[i-1]) - 2.0*h); 151 } 152 153 CHKERRQ(VecRestoreArray(vw,&w)); 154 CHKERRQ(VecRestoreArray(vfw,&fw)); 155 CHKERRQ(DMDAVecRestoreArray(user->da1,vu,&u)); 156 CHKERRQ(DMDAVecRestoreArray(user->da1,vfu,&fu)); 157 CHKERRQ(DMDAVecRestoreArray(user->da1,vlambda,&lambda)); 158 CHKERRQ(DMDAVecRestoreArray(user->da1,vflambda,&flambda)); 159 160 CHKERRQ(DMCompositeGather(user->packer,INSERT_VALUES,FU,vfw,vfu,vflambda)); 161 CHKERRQ(DMCompositeRestoreLocalVectors(user->packer,&vw,&vu,&vlambda)); 162 CHKERRQ(DMCompositeRestoreLocalVectors(user->packer,&vfw,&vfu,&vflambda)); 163 PetscFunctionReturn(0); 164 } 165 166 PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal rnorm,void *dummy) 167 { 168 UserCtx *user = (UserCtx*)dummy; 169 Vec w,u,lambda,U,F; 170 171 PetscFunctionBeginUser; 172 CHKERRQ(SNESGetSolution(snes,&U)); 173 CHKERRQ(DMCompositeGetAccess(user->packer,U,&w,&u,&lambda)); 174 CHKERRQ(VecView(u,user->u_viewer)); 175 CHKERRQ(VecView(lambda,user->lambda_viewer)); 176 CHKERRQ(DMCompositeRestoreAccess(user->packer,U,&w,&u,&lambda)); 177 178 CHKERRQ(SNESGetFunction(snes,&F,0,0)); 179 CHKERRQ(DMCompositeGetAccess(user->packer,F,&w,&u,&lambda)); 180 CHKERRQ(VecView(u,user->fu_viewer)); 181 CHKERRQ(VecView(lambda,user->flambda_viewer)); 182 CHKERRQ(DMCompositeRestoreAccess(user->packer,F,&w,&u,&lambda)); 183 PetscFunctionReturn(0); 184 } 185 186 /*TEST 187 188 test: 189 nsize: 4 190 args: -snes_linesearch_monitor -snes_mf -pc_type none -snes_monitor_short -nox -ksp_monitor_short -snes_converged_reason 191 requires: !single 192 193 TEST*/ 194