Actual source code: ex29.c
1: static char help[] = "Time-Dependent Allan-Cahn example in 2D with Varying Coefficients";
3: /*
4: This example is mainly here to show how to transfer coefficients between subdomains and levels in
5: multigrid and domain decomposition.
6: */
8: #include <petscdm.h>
9: #include <petscdmda.h>
10: #include <petscsnes.h>
11: #include <petscts.h>
13: typedef struct {
14: PetscScalar epsilon;
15: PetscScalar beta;
16: } Coeff;
18: typedef struct {
19: PetscScalar u;
20: } Field;
22: extern PetscErrorCode FormInitialGuess(DM da, void *ctx, Vec X);
23: extern PetscErrorCode FormDiffusionCoefficient(DM da, void *ctx, Vec X);
24: extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo *, PetscReal, Field **, Field **, Field **, void *);
26: /* hooks */
28: static PetscErrorCode CoefficientCoarsenHook(DM dm, DM dmc, void *ctx)
29: {
30: Vec c, cc, ccl;
31: Mat J;
32: Vec vscale;
33: DM cdm, cdmc;
35: PetscFunctionBeginUser;
36: PetscCall(PetscObjectQuery((PetscObject)dm, "coefficientdm", (PetscObject *)&cdm));
38: PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "The coefficient DM needs to be set up!");
40: PetscCall(DMDACreateCompatibleDMDA(dmc, 2, &cdmc));
41: PetscCall(PetscObjectCompose((PetscObject)dmc, "coefficientdm", (PetscObject)cdmc));
43: PetscCall(DMGetNamedGlobalVector(cdm, "coefficient", &c));
44: PetscCall(DMGetNamedGlobalVector(cdmc, "coefficient", &cc));
45: PetscCall(DMGetNamedLocalVector(cdmc, "coefficient", &ccl));
47: PetscCall(DMCreateInterpolation(cdmc, cdm, &J, &vscale));
48: PetscCall(MatRestrict(J, c, cc));
49: PetscCall(VecPointwiseMult(cc, vscale, cc));
51: PetscCall(MatDestroy(&J));
52: PetscCall(VecDestroy(&vscale));
54: PetscCall(DMGlobalToLocalBegin(cdmc, cc, INSERT_VALUES, ccl));
55: PetscCall(DMGlobalToLocalEnd(cdmc, cc, INSERT_VALUES, ccl));
57: PetscCall(DMRestoreNamedGlobalVector(cdm, "coefficient", &c));
58: PetscCall(DMRestoreNamedGlobalVector(cdmc, "coefficient", &cc));
59: PetscCall(DMRestoreNamedLocalVector(cdmc, "coefficient", &ccl));
61: PetscCall(DMCoarsenHookAdd(dmc, CoefficientCoarsenHook, NULL, NULL));
62: PetscCall(DMDestroy(&cdmc));
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: /* This could restrict auxiliary information to the coarse level.
67: */
68: static PetscErrorCode CoefficientSubDomainRestrictHook(DM dm, DM subdm, void *ctx)
69: {
70: Vec c, cc;
71: DM cdm, csubdm;
72: VecScatter *iscat, *oscat, *gscat;
74: PetscFunctionBeginUser;
75: PetscCall(PetscObjectQuery((PetscObject)dm, "coefficientdm", (PetscObject *)&cdm));
77: PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "The coefficient DM needs to be set up!");
79: PetscCall(DMDACreateCompatibleDMDA(subdm, 2, &csubdm));
80: PetscCall(PetscObjectCompose((PetscObject)subdm, "coefficientdm", (PetscObject)csubdm));
82: PetscCall(DMGetNamedGlobalVector(cdm, "coefficient", &c));
83: PetscCall(DMGetNamedLocalVector(csubdm, "coefficient", &cc));
85: PetscCall(DMCreateDomainDecompositionScatters(cdm, 1, &csubdm, &iscat, &oscat, &gscat));
87: PetscCall(VecScatterBegin(*gscat, c, cc, INSERT_VALUES, SCATTER_FORWARD));
88: PetscCall(VecScatterEnd(*gscat, c, cc, INSERT_VALUES, SCATTER_FORWARD));
90: PetscCall(VecScatterDestroy(iscat));
91: PetscCall(VecScatterDestroy(oscat));
92: PetscCall(VecScatterDestroy(gscat));
93: PetscCall(PetscFree(iscat));
94: PetscCall(PetscFree(oscat));
95: PetscCall(PetscFree(gscat));
97: PetscCall(DMRestoreNamedGlobalVector(cdm, "coefficient", &c));
98: PetscCall(DMRestoreNamedLocalVector(csubdm, "coefficient", &cc));
100: PetscCall(DMDestroy(&csubdm));
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: int main(int argc, char **argv)
105: {
106: TS ts;
107: Vec x, c, clocal;
108: DM da, cda;
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));
120: PetscCall(DMDASetFieldName(da, 0, "u"));
121: PetscCall(DMCreateGlobalVector(da, &x));
123: PetscCall(TSSetDM(ts, da));
125: PetscCall(FormInitialGuess(da, NULL, x));
126: PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, PetscReal, void *, void *, void *, void *))FormIFunctionLocal, NULL));
128: /* set up the coefficient */
130: PetscCall(DMDACreateCompatibleDMDA(da, 2, &cda));
131: PetscCall(PetscObjectCompose((PetscObject)da, "coefficientdm", (PetscObject)cda));
133: PetscCall(DMGetNamedGlobalVector(cda, "coefficient", &c));
134: PetscCall(DMGetNamedLocalVector(cda, "coefficient", &clocal));
136: PetscCall(FormDiffusionCoefficient(cda, NULL, c));
138: PetscCall(DMGlobalToLocalBegin(cda, c, INSERT_VALUES, clocal));
139: PetscCall(DMGlobalToLocalEnd(cda, c, INSERT_VALUES, clocal));
141: PetscCall(DMRestoreNamedLocalVector(cda, "coefficient", &clocal));
142: PetscCall(DMRestoreNamedGlobalVector(cda, "coefficient", &c));
144: PetscCall(DMCoarsenHookAdd(da, CoefficientCoarsenHook, NULL, NULL));
145: PetscCall(DMSubDomainHookAdd(da, CoefficientSubDomainRestrictHook, NULL, NULL));
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));
154: PetscCall(TSSolve(ts, x));
156: PetscCall(VecDestroy(&x));
157: PetscCall(TSDestroy(&ts));
158: PetscCall(DMDestroy(&da));
159: PetscCall(DMDestroy(&cda));
161: PetscCall(PetscFinalize());
162: return 0;
163: }
165: /* ------------------------------------------------------------------- */
167: PetscErrorCode FormInitialGuess(DM da, void *ctx, Vec X)
168: {
169: PetscInt i, j, Mx, My, xs, ys, xm, ym;
170: Field **x;
171: PetscReal x0, x1;
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));
176: PetscCall(DMDAVecGetArray(da, X, &x));
177: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
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: }
187: PetscCall(DMDAVecRestoreArray(da, X, &x));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: PetscErrorCode FormDiffusionCoefficient(DM da, void *ctx, Vec X)
192: {
193: PetscInt i, j, Mx, My, xs, ys, xm, ym;
194: Coeff **x;
195: PetscReal x1, x0;
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));
200: /*
201: ierr = VecSetRandom(X,NULL);
202: PetscCall(VecMin(X,NULL,&min));
203: */
205: PetscCall(DMDAVecGetArray(da, X, &x));
206: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
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);
213: x[j][i].epsilon = 0.0;
214: x[j][i].beta = 0.05 + 0.05 * PetscSqrtReal(x0 * x0 + x1 * x1);
215: }
216: }
218: PetscCall(DMDAVecRestoreArray(da, X, &x));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal ptime, Field **x, Field **xt, Field **f, void *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;
232: PetscFunctionBeginUser;
233: PetscCall(PetscObjectQuery((PetscObject)info->da, "coefficientdm", (PetscObject *)&cdm));
234: PetscCall(DMGetNamedLocalVector(cdm, "coefficient", &C));
235: PetscCall(DMDAVecGetArray(cdm, C, &c));
237: hx = 10.0 / ((PetscReal)(info->mx - 1));
238: hy = 10.0 / ((PetscReal)(info->my - 1));
240: dhx = 1. / hx;
241: dhy = 1. / hy;
243: hxdhy = hx / hy;
244: hydhx = hy / hx;
245: scale = hx * hy;
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;
251: u = x[j][i].u;
253: f[j][i].u += scale * (u * u - 1.) * u;
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;
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;
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;
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));
275: PetscCall(DMDAVecRestoreArray(cdm, C, &c));
276: PetscCall(DMRestoreNamedLocalVector(cdm, "coefficient", &C));
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*TEST
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
288: TEST*/