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");PetscCall(ierr); 43 run = options->runType; 44 PetscCall(PetscOptionsEList("-run_type", "The run type", "ex2.c", runTypes, 2, runTypes[options->runType], &run, NULL)); 45 options->runType = (RunType) run; 46 PetscCall(PetscOptionsBool("-use_dual_penalty", "Penalize deviation from both goals", "ex2.c", options->useDualPenalty, &options->useDualPenalty, NULL)); 47 ierr = PetscOptionsEnd();PetscCall(ierr); 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 52 { 53 PetscFunctionBeginUser; 54 PetscCall(DMCreate(comm, dm)); 55 PetscCall(DMSetType(*dm, DMPLEX)); 56 PetscCall(DMSetFromOptions(*dm)); 57 PetscCall(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 PetscCall(DMGetDS(dm, &ds)); 199 PetscCall(PetscDSSetResidual(ds, 0, user->useDualPenalty == PETSC_TRUE ? f0_u_full : f0_u, f1_u)); 200 PetscCall(PetscDSSetResidual(ds, 1, f0_a, NULL)); 201 PetscCall(PetscDSSetResidual(ds, 2, f0_l, f1_l)); 202 PetscCall(PetscDSSetJacobian(ds, 0, 0, user->useDualPenalty == PETSC_TRUE ? g0_uu_full : g0_uu, NULL, NULL, NULL)); 203 PetscCall(PetscDSSetJacobian(ds, 0, 2, NULL, NULL, NULL, g3_ul)); 204 PetscCall(PetscDSSetJacobian(ds, 1, 1, g0_aa, NULL, NULL, NULL)); 205 PetscCall(PetscDSSetJacobian(ds, 1, 2, g0_al, NULL, NULL, NULL)); 206 PetscCall(PetscDSSetJacobian(ds, 2, 1, g0_la, NULL, NULL, NULL)); 207 PetscCall(PetscDSSetJacobian(ds, 2, 0, NULL, NULL, NULL, g3_lu)); 208 209 PetscCall(PetscDSSetExactSolution(ds, 0, quadratic_u_2d, NULL)); 210 PetscCall(PetscDSSetExactSolution(ds, 1, constant_a_2d, NULL)); 211 PetscCall(PetscDSSetExactSolution(ds, 2, zero, NULL)); 212 PetscCall(DMGetLabel(dm, "marker", &label)); 213 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)()) quadratic_u_2d, NULL, user, NULL)); 214 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 1, 0, NULL, (void (*)()) constant_a_2d, NULL, user, NULL)); 215 PetscCall(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 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 230 PetscCall(PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, "potential_", -1, &fe[0])); 231 PetscCall(PetscObjectSetName((PetscObject) fe[0], "potential")); 232 PetscCall(PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, "charge_", -1, &fe[1])); 233 PetscCall(PetscObjectSetName((PetscObject) fe[1], "charge")); 234 PetscCall(PetscFECopyQuadrature(fe[0], fe[1])); 235 PetscCall(PetscFECreateDefault(comm, dim, 1, PETSC_TRUE, "multiplier_", -1, &fe[2])); 236 PetscCall(PetscObjectSetName((PetscObject) fe[2], "multiplier")); 237 PetscCall(PetscFECopyQuadrature(fe[0], fe[2])); 238 /* Set discretization and boundary conditions for each mesh */ 239 for (f = 0; f < 3; ++f) PetscCall(DMSetField(dm, f, NULL, (PetscObject) fe[f])); 240 PetscCall(DMCreateDS(cdm)); 241 PetscCall(SetupProblem(dm, user)); 242 while (cdm) { 243 PetscCall(DMCopyDisc(dm, cdm)); 244 PetscCall(DMGetCoarseDM(cdm, &cdm)); 245 } 246 for (f = 0; f < 3; ++f) PetscCall(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 257 PetscCall(PetscInitialize(&argc, &argv, NULL,help)); 258 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 259 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes)); 260 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm)); 261 PetscCall(SNESSetDM(snes, dm)); 262 PetscCall(SetupDiscretization(dm, &user)); 263 264 PetscCall(DMCreateGlobalVector(dm, &u)); 265 PetscCall(PetscObjectSetName((PetscObject) u, "solution")); 266 PetscCall(VecDuplicate(u, &r)); 267 PetscCall(DMPlexSetSNESLocalFEM(dm,&user,&user,&user)); 268 PetscCall(SNESSetFromOptions(snes)); 269 270 PetscCall(DMSNESCheckFromOptions(snes, u)); 271 if (user.runType == RUN_FULL) { 272 PetscDS ds; 273 PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 274 PetscErrorCode (*initialGuess[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx); 275 PetscReal error; 276 277 PetscCall(DMGetDS(dm, &ds)); 278 PetscCall(PetscDSGetExactSolution(ds, 0, &exactFuncs[0], NULL)); 279 PetscCall(PetscDSGetExactSolution(ds, 1, &exactFuncs[1], NULL)); 280 PetscCall(PetscDSGetExactSolution(ds, 2, &exactFuncs[2], NULL)); 281 initialGuess[0] = zero; 282 initialGuess[1] = zero; 283 initialGuess[2] = zero; 284 PetscCall(DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u)); 285 PetscCall(VecViewFromOptions(u, NULL, "-initial_vec_view")); 286 PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, NULL, u, &error)); 287 if (error < 1.0e-11) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial L_2 Error: < 1.0e-11\n")); 288 else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial L_2 Error: %g\n", error)); 289 PetscCall(SNESSolve(snes, NULL, u)); 290 PetscCall(DMComputeL2Diff(dm, 0.0, exactFuncs, NULL, u, &error)); 291 if (error < 1.0e-11) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final L_2 Error: < 1.0e-11\n")); 292 else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final L_2 Error: %g\n", error)); 293 } 294 PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view")); 295 296 PetscCall(VecDestroy(&u)); 297 PetscCall(VecDestroy(&r)); 298 PetscCall(SNESDestroy(&snes)); 299 PetscCall(DMDestroy(&dm)); 300 PetscCall(PetscFinalize()); 301 return 0; 302 } 303 304 /*TEST 305 306 build: 307 requires: !complex triangle 308 309 test: 310 suffix: 0 311 args: -run_type test -dmsnes_check -potential_petscspace_degree 2 -charge_petscspace_degree 1 -multiplier_petscspace_degree 1 312 313 test: 314 suffix: 1 315 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 316 317 test: 318 suffix: 2 319 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 320 321 TEST*/ 322