1 static char help[] = "One-Shot Multigrid for Parameter Estimation Problem for the Poisson Equation.\n\ 2 Using the Interior Point Method.\n\n\n"; 3 4 /*F 5 We are solving the parameter estimation problem for the Laplacian. We will ask to minimize a Lagrangian 6 function over $y$ and $u$, given by 7 \begin{align} 8 L(u, a, \lambda) = \frac{1}{2} || Qu - d_A ||^2 || Qu - d_B ||^2 + \frac{\beta}{2} || L (a - a_r) ||^2 + \lambda F(u; a) 9 \end{align} 10 where $Q$ is a sampling operator, $L$ is a regularization operator, $F$ defines the PDE. 11 12 Currently, we have perfect information, meaning $Q = I$, and then we need no regularization, $L = I$. We 13 also give the null vector for the reference control $a_r$. Right now $\beta = 1$. 14 15 The PDE will be the Laplace equation with homogeneous boundary conditions 16 \begin{align} 17 -Delta u = a 18 \end{align} 19 20 F*/ 21 22 #include <petsc.h> 23 #include <petscfe.h> 24 25 typedef enum {RUN_FULL, RUN_TEST} RunType; 26 27 typedef struct { 28 RunType runType; /* Whether to run tests, or solve the full problem */ 29 PetscBool useDualPenalty; /* Penalize deviation from both goals */ 30 } AppCtx; 31 32 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 33 { 34 const char *runTypes[2] = {"full", "test"}; 35 PetscInt run; 36 PetscErrorCode ierr; 37 38 PetscFunctionBeginUser; 39 options->runType = RUN_FULL; 40 options->useDualPenalty = PETSC_FALSE; 41 42 ierr = PetscOptionsBegin(comm, "", "Inverse Problem Options", "DMPLEX");CHKERRQ(ierr); 43 run = options->runType; 44 CHKERRQ(PetscOptionsEList("-run_type", "The run type", "ex2.c", runTypes, 2, runTypes[options->runType], &run, NULL)); 45 options->runType = (RunType) run; 46 CHKERRQ(PetscOptionsBool("-use_dual_penalty", "Penalize deviation from both goals", "ex2.c", options->useDualPenalty, &options->useDualPenalty, NULL)); 47 ierr = PetscOptionsEnd();CHKERRQ(ierr); 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 52 { 53 PetscFunctionBeginUser; 54 CHKERRQ(DMCreate(comm, dm)); 55 CHKERRQ(DMSetType(*dm, DMPLEX)); 56 CHKERRQ(DMSetFromOptions(*dm)); 57 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 58 PetscFunctionReturn(0); 59 } 60 61 void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 62 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 63 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 64 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 65 { 66 f0[0] = (u[0] - (x[0]*x[0] + x[1]*x[1])); 67 } 68 void f0_u_full(PetscInt dim, PetscInt Nf, PetscInt NfAux, 69 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 70 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 71 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 72 { 73 f0[0] = (u[0] - (x[0]*x[0] + x[1]*x[1]))*PetscSqr(u[0] - (sin(2.0*PETSC_PI*x[0]) * sin(2.0*PETSC_PI*x[1]))) + 74 PetscSqr(u[0] - (x[0]*x[0] + x[1]*x[1]))*(u[0] - (sin(2.0*PETSC_PI*x[0]) * sin(2.0*PETSC_PI*x[1]))); 75 } 76 void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 77 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 78 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 79 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 80 { 81 PetscInt d; 82 for (d = 0; d < dim; ++d) f1[d] = u_x[dim*2+d]; 83 } 84 void g0_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 85 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 86 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 87 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 88 { 89 g0[0] = 1.0; 90 } 91 void g0_uu_full(PetscInt dim, PetscInt Nf, PetscInt NfAux, 92 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 93 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 94 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 95 { 96 g0[0] = PetscSqr(u[0] - sin(2.0*PETSC_PI*x[0]) * sin(2.0*PETSC_PI*x[1])) 97 + PetscSqr(u[0] - (x[0]*x[0] + x[1]*x[1])) 98 - 2.0*((x[0]*x[0] + x[1]*x[1]) + (sin(2.0*PETSC_PI*x[0]) * sin(2.0*PETSC_PI*x[1])))*u[0] 99 + 4.0*(x[0]*x[0] + x[1]*x[1])*(sin(2.0*PETSC_PI*x[0]) * sin(2.0*PETSC_PI*x[1])); 100 } 101 void g3_ul(PetscInt dim, PetscInt Nf, PetscInt NfAux, 102 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 103 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 104 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 105 { 106 PetscInt d; 107 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 108 } 109 110 void f0_a(PetscInt dim, PetscInt Nf, PetscInt NfAux, 111 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 112 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 113 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 114 { 115 f0[0] = u[1] - 4.0 /* 0.0 */ + u[2]; 116 } 117 void g0_aa(PetscInt dim, PetscInt Nf, PetscInt NfAux, 118 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 119 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 120 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 121 { 122 g0[0] = 1.0; 123 } 124 void g0_al(PetscInt dim, PetscInt Nf, PetscInt NfAux, 125 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 126 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 127 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 128 { 129 g0[0] = 1.0; 130 } 131 132 void f0_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, 133 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 134 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 135 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 136 { 137 f0[0] = u[1]; 138 } 139 void f1_l(PetscInt dim, PetscInt Nf, PetscInt NfAux, 140 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 141 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 142 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 143 { 144 PetscInt d; 145 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 146 } 147 void g0_la(PetscInt dim, PetscInt Nf, PetscInt NfAux, 148 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 149 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 150 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 151 { 152 g0[0] = 1.0; 153 } 154 void g3_lu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 155 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 156 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 157 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 158 { 159 PetscInt d; 160 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 161 } 162 163 /* 164 In 2D for Dirichlet conditions with a variable coefficient, we use exact solution: 165 166 u = x^2 + y^2 167 a = 4 168 d_A = 4 169 d_B = sin(2*pi*x[0]) * sin(2*pi*x[1]) 170 171 so that 172 173 -\Delta u + a = -4 + 4 = 0 174 */ 175 PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 176 { 177 *u = x[0]*x[0] + x[1]*x[1]; 178 return 0; 179 } 180 PetscErrorCode constant_a_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *a, void *ctx) 181 { 182 *a = 4; 183 return 0; 184 } 185 PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *l, void *ctx) 186 { 187 *l = 0.0; 188 return 0; 189 } 190 191 PetscErrorCode SetupProblem(DM dm, AppCtx *user) 192 { 193 PetscDS ds; 194 DMLabel label; 195 const PetscInt id = 1; 196 197 PetscFunctionBeginUser; 198 CHKERRQ(DMGetDS(dm, &ds)); 199 CHKERRQ(PetscDSSetResidual(ds, 0, user->useDualPenalty == PETSC_TRUE ? f0_u_full : f0_u, f1_u)); 200 CHKERRQ(PetscDSSetResidual(ds, 1, f0_a, NULL)); 201 CHKERRQ(PetscDSSetResidual(ds, 2, f0_l, f1_l)); 202 CHKERRQ(PetscDSSetJacobian(ds, 0, 0, user->useDualPenalty == PETSC_TRUE ? g0_uu_full : g0_uu, NULL, NULL, NULL)); 203 CHKERRQ(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, NULL, g3_ul)); 204 CHKERRQ(PetscDSSetJacobian(ds, 1, 1, g0_aa, NULL, NULL, NULL)); 205 CHKERRQ(PetscDSSetJacobian(ds, 1, 2, g0_al, NULL, NULL, NULL)); 206 CHKERRQ(PetscDSSetJacobian(ds, 2, 1, g0_la, NULL, NULL, NULL)); 207 CHKERRQ(PetscDSSetJacobian(ds, 2, 0, NULL, NULL, NULL, g3_lu)); 208 209 CHKERRQ(PetscDSSetExactSolution(ds, 0, quadratic_u_2d, NULL)); 210 CHKERRQ(PetscDSSetExactSolution(ds, 1, constant_a_2d, NULL)); 211 CHKERRQ(PetscDSSetExactSolution(ds, 2, zero, NULL)); 212 CHKERRQ(DMGetLabel(dm, "marker", &label)); 213 CHKERRQ(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)()) quadratic_u_2d, NULL, user, NULL)); 214 CHKERRQ(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 1, 0, NULL, (void (*)()) constant_a_2d, NULL, user, NULL)); 215 CHKERRQ(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 2, 0, NULL, (void (*)()) zero, NULL, user, NULL)); 216 PetscFunctionReturn(0); 217 } 218 219 PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) 220 { 221 DM cdm = dm; 222 const PetscInt dim = 2; 223 PetscFE fe[3]; 224 PetscInt f; 225 MPI_Comm comm; 226 227 PetscFunctionBeginUser; 228 /* Create finite element */ 229 CHKERRQ(PetscObjectGetComm((PetscObject) dm, &comm)); 230 CHKERRQ(PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, "potential_", -1, &fe[0])); 231 CHKERRQ(PetscObjectSetName((PetscObject) fe[0], "potential")); 232 CHKERRQ(PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, "charge_", -1, &fe[1])); 233 CHKERRQ(PetscObjectSetName((PetscObject) fe[1], "charge")); 234 CHKERRQ(PetscFECopyQuadrature(fe[0], fe[1])); 235 CHKERRQ(PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, "multiplier_", -1, &fe[2])); 236 CHKERRQ(PetscObjectSetName((PetscObject) fe[2], "multiplier")); 237 CHKERRQ(PetscFECopyQuadrature(fe[0], fe[2])); 238 /* Set discretization and boundary conditions for each mesh */ 239 for (f = 0; f < 3; ++f) CHKERRQ(DMSetField(dm, f, NULL, (PetscObject) fe[f])); 240 CHKERRQ(DMCreateDS(cdm)); 241 CHKERRQ(SetupProblem(dm, user)); 242 while (cdm) { 243 CHKERRQ(DMCopyDisc(dm, cdm)); 244 CHKERRQ(DMGetCoarseDM(cdm, &cdm)); 245 } 246 for (f = 0; f < 3; ++f) CHKERRQ(PetscFEDestroy(&fe[f])); 247 PetscFunctionReturn(0); 248 } 249 250 int main(int argc, char **argv) 251 { 252 DM dm; 253 SNES snes; 254 Vec u, r; 255 AppCtx user; 256 PetscErrorCode ierr; 257 258 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 259 CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user)); 260 CHKERRQ(SNESCreate(PETSC_COMM_WORLD, &snes)); 261 CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 262 CHKERRQ(SNESSetDM(snes, dm)); 263 CHKERRQ(SetupDiscretization(dm, &user)); 264 265 CHKERRQ(DMCreateGlobalVector(dm, &u)); 266 CHKERRQ(PetscObjectSetName((PetscObject) u, "solution")); 267 CHKERRQ(VecDuplicate(u, &r)); 268 CHKERRQ(DMPlexSetSNESLocalFEM(dm,&user,&user,&user)); 269 CHKERRQ(SNESSetFromOptions(snes)); 270 271 CHKERRQ(DMSNESCheckFromOptions(snes, u)); 272 if (user.runType == RUN_FULL) { 273 PetscDS ds; 274 PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 275 PetscErrorCode (*initialGuess[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx); 276 PetscReal error; 277 278 CHKERRQ(DMGetDS(dm, &ds)); 279 CHKERRQ(PetscDSGetExactSolution(ds, 0, &exactFuncs[0], NULL)); 280 CHKERRQ(PetscDSGetExactSolution(ds, 1, &exactFuncs[1], NULL)); 281 CHKERRQ(PetscDSGetExactSolution(ds, 2, &exactFuncs[2], NULL)); 282 initialGuess[0] = zero; 283 initialGuess[1] = zero; 284 initialGuess[2] = zero; 285 CHKERRQ(DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u)); 286 CHKERRQ(VecViewFromOptions(u, NULL, "-initial_vec_view")); 287 CHKERRQ(DMComputeL2Diff(dm, 0.0, exactFuncs, NULL, u, &error)); 288 if (error < 1.0e-11) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Initial L_2 Error: < 1.0e-11\n")); 289 else CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Initial L_2 Error: %g\n", error)); 290 CHKERRQ(SNESSolve(snes, NULL, u)); 291 CHKERRQ(DMComputeL2Diff(dm, 0.0, exactFuncs, NULL, u, &error)); 292 if (error < 1.0e-11) CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Final L_2 Error: < 1.0e-11\n")); 293 else CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Final L_2 Error: %g\n", error)); 294 } 295 CHKERRQ(VecViewFromOptions(u, NULL, "-sol_vec_view")); 296 297 CHKERRQ(VecDestroy(&u)); 298 CHKERRQ(VecDestroy(&r)); 299 CHKERRQ(SNESDestroy(&snes)); 300 CHKERRQ(DMDestroy(&dm)); 301 ierr = PetscFinalize(); 302 return ierr; 303 } 304 305 /*TEST 306 307 build: 308 requires: !complex triangle 309 310 test: 311 suffix: 0 312 args: -run_type test -dmsnes_check -potential_petscspace_degree 2 -charge_petscspace_degree 1 -multiplier_petscspace_degree 1 313 314 test: 315 suffix: 1 316 args: -potential_petscspace_degree 2 -charge_petscspace_degree 1 -multiplier_petscspace_degree 1 -snes_monitor -snes_converged_reason -pc_type fieldsplit -pc_fieldsplit_0_fields 0,1 -pc_fieldsplit_1_fields 2 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition selfp -fieldsplit_0_pc_type lu -sol_vec_view 317 318 test: 319 suffix: 2 320 args: -potential_petscspace_degree 2 -charge_petscspace_degree 1 -multiplier_petscspace_degree 1 -snes_monitor -snes_converged_reason -snes_fd -pc_type fieldsplit -pc_fieldsplit_0_fields 0,1 -pc_fieldsplit_1_fields 2 -pc_fieldsplit_type schur -pc_fieldsplit_schur_factorization_type full -pc_fieldsplit_schur_precondition selfp -fieldsplit_0_pc_type lu -sol_vec_view 321 322 TEST*/ 323