Actual source code: ex8.c
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";
11: #include <petscdmplex.h>
12: #include <petscsnes.h>
13: #include <petscds.h>
15: /* The pointwise evaluation routine for the objective function */
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]);
21: obj[0] = PetscSqrtReal(PetscAbsScalar(1 + ux2 + uy2));
22: }
24: /* The pointwise evaluation routine for the gradient wrt the gradients of the FEM basis */
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));
31: f[0] = v * u_x[0];
32: f[1] = v * u_x[1];
33: }
35: /* The pointwise evaluation routine for the hessian wrt the gradients of the FEM basis */
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);
44: jac[0] = v1 - v2 * ux2;
45: jac[1] = -v2 * uxy;
46: jac[2] = -v2 * uxy;
47: jac[3] = v1 - v2 * uy2;
48: }
50: /* The boundary conditions we impose */
51: static PetscErrorCode sins_2d(PetscInt dim, PetscReal time, const PetscReal xx[], PetscInt Nc, PetscScalar *u, void *ctx)
52: {
53: PetscFunctionBeginUser;
54: const PetscReal pi2 = PETSC_PI * 2;
55: const PetscReal x = xx[0];
56: const PetscReal y = xx[1];
58: *u = (y - 0.5) * PetscSinReal(pi2 * x) + (x - 0.5) * PetscSinReal(pi2 * y);
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
63: {
64: DM plex;
65: DMLabel label;
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: }
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: }
86: static PetscErrorCode SetupProblem(DM dm)
87: {
88: PetscDS ds;
89: DMLabel label;
90: PetscInt dim;
91: const PetscInt id = 1;
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, (void (*)(void))sins_2d, NULL, NULL, NULL));
102: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_TRUE, NULL));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: static PetscErrorCode SetupDiscretization(DM dm)
107: {
108: DM plex, cdm = dm;
109: PetscFE fe;
110: PetscBool simplex;
111: PetscInt dim;
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;
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: }
135: int main(int argc, char **argv)
136: {
137: DM dm; /* Problem specification */
138: SNES snes; /* nonlinear solver */
139: Vec u; /* solution vector */
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));
147: PetscCall(SetupDiscretization(dm));
149: PetscCall(SNESSetFromOptions(snes));
151: PetscCall(DMCreateGlobalVector(dm, &u));
152: PetscCall(SNESSolve(snes, NULL, u));
154: PetscCall(VecDestroy(&u));
155: PetscCall(SNESDestroy(&snes));
156: PetscCall(DMDestroy(&dm));
157: PetscCall(PetscFinalize());
158: return 0;
159: }
161: /*TEST
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
170: TEST*/