57 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)  in CreateMesh()  argument
64   user->use_riesz = PETSC_TRUE;  in CreateMesh()
67 …esz", "Use the Riesz map to achieve mesh independence", "ex3.c", user->use_riesz, &user->use_riesz…  in CreateMesh()
177 PetscErrorCode CreateCtx(DM dm, AppCtx *user)  in CreateCtx()  argument
199   PetscCall(DMCreateGlobalVector(dm, &user->data));  in CreateCtx()
205   PetscCall(DMProjectFunction(dm, 0.0, wtf, NULL, INSERT_VALUES, user->data));  in CreateCtx()
216   PetscCall(DMCreateMatrix(dm_mass, &user->mass));  in CreateCtx()
217   PetscCall(DMPlexSNESComputeJacobianFEM(dm_mass, user->data, user->mass, user->mass, NULL));  in CreateCtx()
218   PetscCall(MatSetOption(user->mass, MAT_SYMMETRIC, PETSC_TRUE));  in CreateCtx()
229   PetscCall(DMCreateMatrix(dm_laplace, &user->laplace));  in CreateCtx()
230 …PetscCall(DMPlexSNESComputeJacobianFEM(dm_laplace, user->data, user->laplace, user->laplace, NULL)…  in CreateCtx()
237   user->num_bc_dofs = 0;  in CreateCtx()
240     user->num_bc_dofs += dof;  in CreateCtx()
242   PetscCall(PetscMalloc1(user->num_bc_dofs, &user->bc_indices));  in CreateCtx()
246     for (d = 0; d < dof; ++d) user->bc_indices[k++] = off + d;  in CreateCtx()
255   PetscCall(MatZeroRows(user->laplace, user->num_bc_dofs, user->bc_indices, 1.0, NULL, NULL));  in CreateCtx()
256   PetscCall(PetscCalloc1(user->num_bc_dofs, &user->bc_values));  in CreateCtx()
259   PetscCall(KSPCreate(PETSC_COMM_WORLD, &user->ksp_laplace));  in CreateCtx()
260   PetscCall(KSPSetOperators(user->ksp_laplace, user->laplace, user->laplace));  in CreateCtx()
261   PetscCall(KSPSetOptionsPrefix(user->ksp_laplace, "laplace_"));  in CreateCtx()
262   PetscCall(KSPSetFromOptions(user->ksp_laplace));  in CreateCtx()
265   user->dm = dm;  in CreateCtx()
266   PetscCall(VecDuplicate(user->data, &user->state));  in CreateCtx()
267   PetscCall(VecDuplicate(user->data, &user->adjoint));  in CreateCtx()
268   PetscCall(VecDuplicate(user->data, &user->tmp1));  in CreateCtx()
269   PetscCall(VecDuplicate(user->data, &user->tmp2));  in CreateCtx()
275 PetscErrorCode DestroyCtx(AppCtx *user)  in DestroyCtx()  argument
278   PetscCall(MatDestroy(&user->mass));  in DestroyCtx()
279   PetscCall(MatDestroy(&user->laplace));  in DestroyCtx()
280   PetscCall(KSPDestroy(&user->ksp_laplace));  in DestroyCtx()
281   PetscCall(VecDestroy(&user->data));  in DestroyCtx()
282   PetscCall(VecDestroy(&user->state));  in DestroyCtx()
283   PetscCall(VecDestroy(&user->adjoint));  in DestroyCtx()
284   PetscCall(VecDestroy(&user->tmp1));  in DestroyCtx()
285   PetscCall(VecDestroy(&user->tmp2));  in DestroyCtx()
286   PetscCall(PetscFree(user->bc_indices));  in DestroyCtx()
287   PetscCall(PetscFree(user->bc_values));  in DestroyCtx()
293   AppCtx         *user  = (AppCtx *)userv;  in ReducedFunctionGradient()  local
298   PetscCall(MatMult(user->mass, u, user->tmp1));  in ReducedFunctionGradient()
299   PetscCall(VecDot(u, user->tmp1, &inner)); /* regularisation contribution to */  in ReducedFunctionGradient()
303   PetscCall(VecAXPY(g, alpha, user->tmp1)); /* regularisation contribution to the gradient */  in ReducedFunctionGradient()
306 …PetscCall(VecSetValues(user->tmp1, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VA…  in ReducedFunctionGradient()
307   PetscCall(VecAssemblyBegin(user->tmp1));  in ReducedFunctionGradient()
308   PetscCall(VecAssemblyEnd(user->tmp1));  in ReducedFunctionGradient()
309   PetscCall(KSPSolve(user->ksp_laplace, user->tmp1, user->state)); /* forward solve */  in ReducedFunctionGradient()
312   PetscCall(VecCopy(user->state, user->tmp1));  in ReducedFunctionGradient()
313   PetscCall(VecAXPY(user->tmp1, -1.0, user->data));  in ReducedFunctionGradient()
314   PetscCall(MatMult(user->mass, user->tmp1, user->tmp2));  in ReducedFunctionGradient()
315   PetscCall(VecDot(user->tmp1, user->tmp2, &inner)); /* misfit contribution to */  in ReducedFunctionGradient()
318 …PetscCall(VecSetValues(user->tmp2, user->num_bc_dofs, user->bc_indices, user->bc_values, INSERT_VA…  in ReducedFunctionGradient()
319   PetscCall(VecAssemblyBegin(user->tmp2));  in ReducedFunctionGradient()
320   PetscCall(VecAssemblyEnd(user->tmp2));  in ReducedFunctionGradient()
321   PetscCall(KSPSolve(user->ksp_laplace, user->tmp2, user->adjoint)); /* adjoint solve */  in ReducedFunctionGradient()
324   PetscCall(MatMult(user->mass, user->adjoint, user->tmp1));  in ReducedFunctionGradient()
325   PetscCall(VecAXPY(g, 1.0, user->tmp1)); /* adjoint contribution to the gradient */  in ReducedFunctionGradient()
334   AppCtx user;  in main()  local
338   PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));  in main()
339   PetscCall(CreateCtx(dm, &user));  in main()
350   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, ReducedFunctionGradient, &user));  in main()
355   if (user.use_riesz) {  in main()
356     PetscCall(TaoLMVMSetH0(tao, user.mass)); /* crucial for mesh independence */  in main()
357     PetscCall(TaoSetGradientNorm(tao, user.mass));  in main()
370   PetscCall(DestroyCtx(&user));  in main()