1 static char help[] = "Time-Dependent Allan-Cahn example in 2D with Varying Coefficients";
2
3 /*
4 This example is mainly here to show how to transfer coefficients between subdomains and levels in
5 multigrid and domain decomposition.
6 */
7
8 #include <petscdm.h>
9 #include <petscdmda.h>
10 #include <petscsnes.h>
11 #include <petscts.h>
12
13 typedef struct {
14 PetscScalar epsilon;
15 PetscScalar beta;
16 } Coeff;
17
18 typedef struct {
19 PetscScalar u;
20 } Field;
21
22 extern PetscErrorCode FormInitialGuess(DM da, PetscCtx ctx, Vec X);
23 extern PetscErrorCode FormDiffusionCoefficient(DM da, PetscCtx ctx, Vec X);
24 extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo *, PetscReal, Field **, Field **, Field **, void *);
25
26 /* hooks */
27
CoefficientCoarsenHook(DM dm,DM dmc,PetscCtx ctx)28 static PetscErrorCode CoefficientCoarsenHook(DM dm, DM dmc, PetscCtx ctx)
29 {
30 Vec c, cc, ccl;
31 Mat J;
32 Vec vscale;
33 DM cdm, cdmc;
34
35 PetscFunctionBeginUser;
36 PetscCall(PetscObjectQuery((PetscObject)dm, "coefficientdm", (PetscObject *)&cdm));
37
38 PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "The coefficient DM needs to be set up!");
39
40 PetscCall(DMDACreateCompatibleDMDA(dmc, 2, &cdmc));
41 PetscCall(PetscObjectCompose((PetscObject)dmc, "coefficientdm", (PetscObject)cdmc));
42
43 PetscCall(DMGetNamedGlobalVector(cdm, "coefficient", &c));
44 PetscCall(DMGetNamedGlobalVector(cdmc, "coefficient", &cc));
45 PetscCall(DMGetNamedLocalVector(cdmc, "coefficient", &ccl));
46
47 PetscCall(DMCreateInterpolation(cdmc, cdm, &J, &vscale));
48 PetscCall(MatRestrict(J, c, cc));
49 PetscCall(VecPointwiseMult(cc, vscale, cc));
50
51 PetscCall(MatDestroy(&J));
52 PetscCall(VecDestroy(&vscale));
53
54 PetscCall(DMGlobalToLocalBegin(cdmc, cc, INSERT_VALUES, ccl));
55 PetscCall(DMGlobalToLocalEnd(cdmc, cc, INSERT_VALUES, ccl));
56
57 PetscCall(DMRestoreNamedGlobalVector(cdm, "coefficient", &c));
58 PetscCall(DMRestoreNamedGlobalVector(cdmc, "coefficient", &cc));
59 PetscCall(DMRestoreNamedLocalVector(cdmc, "coefficient", &ccl));
60
61 PetscCall(DMCoarsenHookAdd(dmc, CoefficientCoarsenHook, NULL, NULL));
62 PetscCall(DMDestroy(&cdmc));
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
66 /* This could restrict auxiliary information to the coarse level.
67 */
CoefficientSubDomainRestrictHook(DM dm,DM subdm,PetscCtx ctx)68 static PetscErrorCode CoefficientSubDomainRestrictHook(DM dm, DM subdm, PetscCtx ctx)
69 {
70 Vec c, cc;
71 DM cdm, csubdm;
72 VecScatter *iscat, *oscat, *gscat;
73
74 PetscFunctionBeginUser;
75 PetscCall(PetscObjectQuery((PetscObject)dm, "coefficientdm", (PetscObject *)&cdm));
76
77 PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "The coefficient DM needs to be set up!");
78
79 PetscCall(DMDACreateCompatibleDMDA(subdm, 2, &csubdm));
80 PetscCall(PetscObjectCompose((PetscObject)subdm, "coefficientdm", (PetscObject)csubdm));
81
82 PetscCall(DMGetNamedGlobalVector(cdm, "coefficient", &c));
83 PetscCall(DMGetNamedLocalVector(csubdm, "coefficient", &cc));
84
85 PetscCall(DMCreateDomainDecompositionScatters(cdm, 1, &csubdm, &iscat, &oscat, &gscat));
86
87 PetscCall(VecScatterBegin(*gscat, c, cc, INSERT_VALUES, SCATTER_FORWARD));
88 PetscCall(VecScatterEnd(*gscat, c, cc, INSERT_VALUES, SCATTER_FORWARD));
89
90 PetscCall(VecScatterDestroy(iscat));
91 PetscCall(VecScatterDestroy(oscat));
92 PetscCall(VecScatterDestroy(gscat));
93 PetscCall(PetscFree(iscat));
94 PetscCall(PetscFree(oscat));
95 PetscCall(PetscFree(gscat));
96
97 PetscCall(DMRestoreNamedGlobalVector(cdm, "coefficient", &c));
98 PetscCall(DMRestoreNamedLocalVector(csubdm, "coefficient", &cc));
99
100 PetscCall(DMDestroy(&csubdm));
101 PetscFunctionReturn(PETSC_SUCCESS);
102 }
103
main(int argc,char ** argv)104 int main(int argc, char **argv)
105 {
106 TS ts;
107 Vec x, c, clocal;
108 DM da, cda;
109
110 PetscFunctionBeginUser;
111 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
112 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
113 PetscCall(TSSetType(ts, TSARKIMEX));
114 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
115 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
116 PetscCall(DMSetFromOptions(da));
117 PetscCall(DMSetUp(da));
118 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
119
120 PetscCall(DMDASetFieldName(da, 0, "u"));
121 PetscCall(DMCreateGlobalVector(da, &x));
122
123 PetscCall(TSSetDM(ts, da));
124
125 PetscCall(FormInitialGuess(da, NULL, x));
126 PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, PetscReal, void *, void *, void *, void *))FormIFunctionLocal, NULL));
127
128 /* set up the coefficient */
129
130 PetscCall(DMDACreateCompatibleDMDA(da, 2, &cda));
131 PetscCall(PetscObjectCompose((PetscObject)da, "coefficientdm", (PetscObject)cda));
132
133 PetscCall(DMGetNamedGlobalVector(cda, "coefficient", &c));
134 PetscCall(DMGetNamedLocalVector(cda, "coefficient", &clocal));
135
136 PetscCall(FormDiffusionCoefficient(cda, NULL, c));
137
138 PetscCall(DMGlobalToLocalBegin(cda, c, INSERT_VALUES, clocal));
139 PetscCall(DMGlobalToLocalEnd(cda, c, INSERT_VALUES, clocal));
140
141 PetscCall(DMRestoreNamedLocalVector(cda, "coefficient", &clocal));
142 PetscCall(DMRestoreNamedGlobalVector(cda, "coefficient", &c));
143
144 PetscCall(DMCoarsenHookAdd(da, CoefficientCoarsenHook, NULL, NULL));
145 PetscCall(DMSubDomainHookAdd(da, CoefficientSubDomainRestrictHook, NULL, NULL));
146
147 PetscCall(TSSetMaxSteps(ts, 10000));
148 PetscCall(TSSetMaxTime(ts, 10000.0));
149 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
150 PetscCall(TSSetTimeStep(ts, 0.05));
151 PetscCall(TSSetSolution(ts, x));
152 PetscCall(TSSetFromOptions(ts));
153
154 PetscCall(TSSolve(ts, x));
155
156 PetscCall(VecDestroy(&x));
157 PetscCall(TSDestroy(&ts));
158 PetscCall(DMDestroy(&da));
159 PetscCall(DMDestroy(&cda));
160
161 PetscCall(PetscFinalize());
162 return 0;
163 }
164
165 /* ------------------------------------------------------------------- */
166
FormInitialGuess(DM da,PetscCtx ctx,Vec X)167 PetscErrorCode FormInitialGuess(DM da, PetscCtx ctx, Vec X)
168 {
169 PetscInt i, j, Mx, My, xs, ys, xm, ym;
170 Field **x;
171 PetscReal x0, x1;
172
173 PetscFunctionBeginUser;
174 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
175
176 PetscCall(DMDAVecGetArray(da, X, &x));
177 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
178
179 for (j = ys; j < ys + ym; j++) {
180 for (i = xs; i < xs + xm; i++) {
181 x0 = 10.0 * (i - 0.5 * (Mx - 1)) / (Mx - 1);
182 x1 = 10.0 * (j - 0.5 * (Mx - 1)) / (My - 1);
183 x[j][i].u = PetscCosReal(2.0 * PetscSqrtReal(x1 * x1 + x0 * x0));
184 }
185 }
186
187 PetscCall(DMDAVecRestoreArray(da, X, &x));
188 PetscFunctionReturn(PETSC_SUCCESS);
189 }
190
FormDiffusionCoefficient(DM da,PetscCtx ctx,Vec X)191 PetscErrorCode FormDiffusionCoefficient(DM da, PetscCtx ctx, Vec X)
192 {
193 PetscInt i, j, Mx, My, xs, ys, xm, ym;
194 Coeff **x;
195 PetscReal x1, x0;
196
197 PetscFunctionBeginUser;
198 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
199
200 /*
201 ierr = VecSetRandom(X,NULL);
202 PetscCall(VecMin(X,NULL,&min));
203 */
204
205 PetscCall(DMDAVecGetArray(da, X, &x));
206 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
207
208 for (j = ys; j < ys + ym; j++) {
209 for (i = xs; i < xs + xm; i++) {
210 x0 = 10.0 * (i - 0.5 * (Mx - 1)) / (Mx - 1);
211 x1 = 10.0 * (j - 0.5 * (My - 1)) / (My - 1);
212
213 x[j][i].epsilon = 0.0;
214 x[j][i].beta = 0.05 + 0.05 * PetscSqrtReal(x0 * x0 + x1 * x1);
215 }
216 }
217
218 PetscCall(DMDAVecRestoreArray(da, X, &x));
219 PetscFunctionReturn(PETSC_SUCCESS);
220 }
221
FormIFunctionLocal(DMDALocalInfo * info,PetscReal ptime,Field ** x,Field ** xt,Field ** f,PetscCtx ctx)222 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal ptime, Field **x, Field **xt, Field **f, PetscCtx ctx)
223 {
224 PetscInt i, j;
225 PetscReal hx, hy, dhx, dhy, hxdhy, hydhx, scale;
226 PetscScalar u, uxx, uyy;
227 PetscScalar ux, uy, bx, by;
228 Vec C;
229 Coeff **c;
230 DM cdm;
231
232 PetscFunctionBeginUser;
233 PetscCall(PetscObjectQuery((PetscObject)info->da, "coefficientdm", (PetscObject *)&cdm));
234 PetscCall(DMGetNamedLocalVector(cdm, "coefficient", &C));
235 PetscCall(DMDAVecGetArray(cdm, C, &c));
236
237 hx = 10.0 / ((PetscReal)(info->mx - 1));
238 hy = 10.0 / ((PetscReal)(info->my - 1));
239
240 dhx = 1. / hx;
241 dhy = 1. / hy;
242
243 hxdhy = hx / hy;
244 hydhx = hy / hx;
245 scale = hx * hy;
246
247 for (j = info->ys; j < info->ys + info->ym; j++) {
248 for (i = info->xs; i < info->xs + info->xm; i++) {
249 f[j][i].u = xt[j][i].u * scale;
250
251 u = x[j][i].u;
252
253 f[j][i].u += scale * (u * u - 1.) * u;
254
255 if (i == 0) f[j][i].u += (x[j][i].u - x[j][i + 1].u) * dhx;
256 else if (i == info->mx - 1) f[j][i].u += (x[j][i].u - x[j][i - 1].u) * dhx;
257 else if (j == 0) f[j][i].u += (x[j][i].u - x[j + 1][i].u) * dhy;
258 else if (j == info->my - 1) f[j][i].u += (x[j][i].u - x[j - 1][i].u) * dhy;
259 else {
260 uyy = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;
261 uxx = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;
262
263 bx = 0.5 * (c[j][i + 1].beta - c[j][i - 1].beta) * dhx;
264 by = 0.5 * (c[j + 1][i].beta - c[j - 1][i].beta) * dhy;
265
266 ux = 0.5 * (x[j][i + 1].u - x[j][i - 1].u) * dhx;
267 uy = 0.5 * (x[j + 1][i].u - x[j - 1][i].u) * dhy;
268
269 f[j][i].u += c[j][i].beta * (uxx + uyy) + scale * (bx * ux + by * uy);
270 }
271 }
272 }
273 PetscCall(PetscLogFlops(11. * info->ym * info->xm));
274
275 PetscCall(DMDAVecRestoreArray(cdm, C, &c));
276 PetscCall(DMRestoreNamedLocalVector(cdm, "coefficient", &C));
277 PetscFunctionReturn(PETSC_SUCCESS);
278 }
279
280 /*TEST
281
282 test:
283 args: -da_refine 4 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 -ts_type arkimex -ts_monitor -snes_monitor -snes_type ngmres -npc_snes_type nasm -npc_snes_nasm_type restrict -da_overlap 4
284 nsize: 16
285 requires: !single
286 output_file: output/ex29.out
287
288 TEST*/
289