1 static const char help[] = "\
2 Minimum surface area problem in 2D using DMPLEX.\n\
3 It solves the unconstrained minimization problem \n\
4 \n\
5 argmin \\int_\\Omega (1 + ||u||^2), u = g on \\partial\\Omega.\n\
6 \n\
7 This example shows how to setup an optimization problem using DMPLEX FEM routines.\n\
8 It supports nonlinear domain decomposition and multilevel solvers.\n\
9 \n\n";
10
11 #include <petscdmplex.h>
12 #include <petscsnes.h>
13 #include <petscds.h>
14
15 /* The pointwise evaluation routine for the objective function */
objective_2d(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar obj[])16 static void objective_2d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
17 {
18 const PetscScalar ux2 = PetscSqr(u_x[0]);
19 const PetscScalar uy2 = PetscSqr(u_x[1]);
20
21 obj[0] = PetscSqrtReal(PetscAbsScalar(1 + ux2 + uy2));
22 }
23
24 /* The pointwise evaluation routine for the gradient wrt the gradients of the FEM basis */
gradient_1_2d(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f[])25 static void gradient_1_2d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[])
26 {
27 const PetscScalar ux2 = PetscSqr(u_x[0]);
28 const PetscScalar uy2 = PetscSqr(u_x[1]);
29 const PetscScalar v = 1. / PetscSqrtReal(PetscAbsScalar(1 + ux2 + uy2));
30
31 f[0] = v * u_x[0];
32 f[1] = v * u_x[1];
33 }
34
35 /* The pointwise evaluation routine for the hessian wrt the gradients of the FEM basis */
hessian_11_2d(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar jac[])36 static void hessian_11_2d(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar jac[])
37 {
38 const PetscScalar ux2 = PetscSqr(u_x[0]);
39 const PetscScalar uy2 = PetscSqr(u_x[1]);
40 const PetscScalar uxy = u_x[0] * u_x[1];
41 const PetscScalar v1 = 1. / PetscSqrtReal(PetscAbsScalar(1 + ux2 + uy2));
42 const PetscScalar v2 = v1 / (1 + ux2 + uy2);
43
44 jac[0] = v1 - v2 * ux2;
45 jac[1] = -v2 * uxy;
46 jac[2] = -v2 * uxy;
47 jac[3] = v1 - v2 * uy2;
48 }
49
50 /* The boundary conditions we impose */
sins_2d(PetscInt dim,PetscReal time,const PetscReal xx[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)51 static PetscErrorCode sins_2d(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
52 {
53 PetscFunctionBeginUser;
54 const PetscReal pi2 = PETSC_PI * 2;
55 const PetscReal x = xx[0];
56 const PetscReal y = xx[1];
57
58 *u = (y - 0.5) * PetscSinReal(pi2 * x) + (x - 0.5) * PetscSinReal(pi2 * y);
59 PetscFunctionReturn(PETSC_SUCCESS);
60 }
61
CreateBCLabel(DM dm,const char name[])62 static PetscErrorCode CreateBCLabel(DM dm, const char name[])
63 {
64 DM plex;
65 DMLabel label;
66
67 PetscFunctionBeginUser;
68 PetscCall(DMCreateLabel(dm, name));
69 PetscCall(DMGetLabel(dm, name, &label));
70 PetscCall(DMConvert(dm, DMPLEX, &plex));
71 PetscCall(DMPlexMarkBoundaryFaces(plex, 1, label));
72 PetscCall(DMDestroy(&plex));
73 PetscFunctionReturn(PETSC_SUCCESS);
74 }
75
CreateMesh(MPI_Comm comm,DM * dm)76 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm)
77 {
78 PetscFunctionBeginUser;
79 PetscCall(DMCreate(comm, dm));
80 PetscCall(DMSetType(*dm, DMPLEX));
81 PetscCall(DMSetFromOptions(*dm));
82 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85
SetupProblem(DM dm)86 static PetscErrorCode SetupProblem(DM dm)
87 {
88 PetscDS ds;
89 DMLabel label;
90 PetscInt dim;
91 const PetscInt id = 1;
92
93 PetscFunctionBeginUser;
94 PetscCall(DMGetDS(dm, &ds));
95 PetscCall(DMGetDimension(dm, &dim));
96 PetscCheck(dim == 2, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only for 2-D");
97 PetscCall(PetscDSSetObjective(ds, 0, objective_2d));
98 PetscCall(PetscDSSetResidual(ds, 0, NULL, gradient_1_2d));
99 PetscCall(PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, hessian_11_2d));
100 PetscCall(DMGetLabel(dm, "marker", &label));
101 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "data", label, 1, &id, 0, 0, NULL, (PetscVoidFn *)sins_2d, NULL, NULL, NULL));
102 PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_TRUE, NULL));
103 PetscFunctionReturn(PETSC_SUCCESS);
104 }
105
SetupDiscretization(DM dm)106 static PetscErrorCode SetupDiscretization(DM dm)
107 {
108 DM plex, cdm = dm;
109 PetscFE fe;
110 PetscBool simplex;
111 PetscInt dim;
112
113 PetscFunctionBeginUser;
114 PetscCall(DMGetDimension(dm, &dim));
115 PetscCall(DMConvert(dm, DMPLEX, &plex));
116 PetscCall(DMPlexIsSimplex(plex, &simplex));
117 PetscCall(DMDestroy(&plex));
118 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, NULL, -1, &fe));
119 PetscCall(PetscObjectSetName((PetscObject)fe, "potential"));
120 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe));
121 PetscCall(DMCreateDS(dm));
122 PetscCall(SetupProblem(dm));
123 while (cdm) {
124 PetscBool hasLabel;
125
126 PetscCall(DMHasLabel(cdm, "marker", &hasLabel));
127 if (!hasLabel) PetscCall(CreateBCLabel(cdm, "marker"));
128 PetscCall(DMCopyDisc(dm, cdm));
129 PetscCall(DMGetCoarseDM(cdm, &cdm));
130 }
131 PetscCall(PetscFEDestroy(&fe));
132 PetscFunctionReturn(PETSC_SUCCESS);
133 }
134
main(int argc,char ** argv)135 int main(int argc, char **argv)
136 {
137 DM dm; /* Problem specification */
138 SNES snes; /* nonlinear solver */
139 Vec u; /* solution vector */
140
141 PetscFunctionBeginUser;
142 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
143 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
144 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm));
145 PetscCall(SNESSetDM(snes, dm));
146
147 PetscCall(SetupDiscretization(dm));
148
149 PetscCall(SNESSetFromOptions(snes));
150
151 PetscCall(DMCreateGlobalVector(dm, &u));
152 PetscCall(SNESSolve(snes, NULL, u));
153
154 PetscCall(VecDestroy(&u));
155 PetscCall(SNESDestroy(&snes));
156 PetscCall(DMDestroy(&dm));
157 PetscCall(PetscFinalize());
158 return 0;
159 }
160
161 /*TEST
162
163 test:
164 requires: !single
165 suffix: qn_nasm
166 filter: sed -e "s/CONVERGED_FNORM_ABS/CONVERGED_FNORM_RELATIVE/g"
167 nsize: 4
168 args: -petscspace_degree 1 -dm_refine 2 -snes_type qn -snes_npc_side {{left right}separate output} -npc_snes_type nasm -npc_snes_nasm_type restrict -npc_sub_snes_linesearch_order 1 -npc_sub_snes_linesearch_type bt -dm_plex_dd_overlap 1 -snes_linesearch_type bt -snes_linesearch_order 1 -npc_sub_pc_type lu -npc_sub_ksp_type preonly -snes_converged_reason -snes_monitor_short -petscpartitioner_type simple -npc_sub_snes_max_it 1 -dm_plex_simplex 0 -snes_rtol 1.e-6
169
170 TEST*/
171