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;
36: PetscObjectQuery((PetscObject)dm, "coefficientdm", (PetscObject *)&cdm);
40: DMDACreateCompatibleDMDA(dmc, 2, &cdmc);
41: PetscObjectCompose((PetscObject)dmc, "coefficientdm", (PetscObject)cdmc);
43: DMGetNamedGlobalVector(cdm, "coefficient", &c);
44: DMGetNamedGlobalVector(cdmc, "coefficient", &cc);
45: DMGetNamedLocalVector(cdmc, "coefficient", &ccl);
47: DMCreateInterpolation(cdmc, cdm, &J, &vscale);
48: MatRestrict(J, c, cc);
49: VecPointwiseMult(cc, vscale, cc);
51: MatDestroy(&J);
52: VecDestroy(&vscale);
54: DMGlobalToLocalBegin(cdmc, cc, INSERT_VALUES, ccl);
55: DMGlobalToLocalEnd(cdmc, cc, INSERT_VALUES, ccl);
57: DMRestoreNamedGlobalVector(cdm, "coefficient", &c);
58: DMRestoreNamedGlobalVector(cdmc, "coefficient", &cc);
59: DMRestoreNamedLocalVector(cdmc, "coefficient", &ccl);
61: DMCoarsenHookAdd(dmc, CoefficientCoarsenHook, NULL, NULL);
62: DMDestroy(&cdmc);
63: return 0;
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;
75: PetscObjectQuery((PetscObject)dm, "coefficientdm", (PetscObject *)&cdm);
79: DMDACreateCompatibleDMDA(subdm, 2, &csubdm);
80: PetscObjectCompose((PetscObject)subdm, "coefficientdm", (PetscObject)csubdm);
82: DMGetNamedGlobalVector(cdm, "coefficient", &c);
83: DMGetNamedLocalVector(csubdm, "coefficient", &cc);
85: DMCreateDomainDecompositionScatters(cdm, 1, &csubdm, &iscat, &oscat, &gscat);
87: VecScatterBegin(*gscat, c, cc, INSERT_VALUES, SCATTER_FORWARD);
88: VecScatterEnd(*gscat, c, cc, INSERT_VALUES, SCATTER_FORWARD);
90: VecScatterDestroy(iscat);
91: VecScatterDestroy(oscat);
92: VecScatterDestroy(gscat);
93: PetscFree(iscat);
94: PetscFree(oscat);
95: PetscFree(gscat);
97: DMRestoreNamedGlobalVector(cdm, "coefficient", &c);
98: DMRestoreNamedLocalVector(csubdm, "coefficient", &cc);
100: DMDestroy(&csubdm);
101: return 0;
102: }
104: int main(int argc, char **argv)
106: {
107: TS ts;
108: Vec x, c, clocal;
109: DM da, cda;
112: PetscInitialize(&argc, &argv, (char *)0, help);
113: TSCreate(PETSC_COMM_WORLD, &ts);
114: TSSetType(ts, TSARKIMEX);
115: TSSetProblemType(ts, TS_NONLINEAR);
116: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da);
117: DMSetFromOptions(da);
118: DMSetUp(da);
119: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
121: DMDASetFieldName(da, 0, "u");
122: DMCreateGlobalVector(da, &x);
124: TSSetDM(ts, da);
126: FormInitialGuess(da, NULL, x);
127: DMDATSSetIFunctionLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, PetscReal, void *, void *, void *, void *))FormIFunctionLocal, NULL);
129: /* set up the coefficient */
131: DMDACreateCompatibleDMDA(da, 2, &cda);
132: PetscObjectCompose((PetscObject)da, "coefficientdm", (PetscObject)cda);
134: DMGetNamedGlobalVector(cda, "coefficient", &c);
135: DMGetNamedLocalVector(cda, "coefficient", &clocal);
137: FormDiffusionCoefficient(cda, NULL, c);
139: DMGlobalToLocalBegin(cda, c, INSERT_VALUES, clocal);
140: DMGlobalToLocalEnd(cda, c, INSERT_VALUES, clocal);
142: DMRestoreNamedLocalVector(cda, "coefficient", &clocal);
143: DMRestoreNamedGlobalVector(cda, "coefficient", &c);
145: DMCoarsenHookAdd(da, CoefficientCoarsenHook, NULL, NULL);
146: DMSubDomainHookAdd(da, CoefficientSubDomainRestrictHook, NULL, NULL);
148: TSSetMaxSteps(ts, 10000);
149: TSSetMaxTime(ts, 10000.0);
150: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
151: TSSetTimeStep(ts, 0.05);
152: TSSetSolution(ts, x);
153: TSSetFromOptions(ts);
155: TSSolve(ts, x);
157: VecDestroy(&x);
158: TSDestroy(&ts);
159: DMDestroy(&da);
160: DMDestroy(&cda);
162: PetscFinalize();
163: return 0;
164: }
166: /* ------------------------------------------------------------------- */
168: PetscErrorCode FormInitialGuess(DM da, void *ctx, Vec X)
169: {
170: PetscInt i, j, Mx, My, xs, ys, xm, ym;
171: Field **x;
172: PetscReal x0, x1;
175: 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);
177: DMDAVecGetArray(da, X, &x);
178: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
180: for (j = ys; j < ys + ym; j++) {
181: for (i = xs; i < xs + xm; i++) {
182: x0 = 10.0 * (i - 0.5 * (Mx - 1)) / (Mx - 1);
183: x1 = 10.0 * (j - 0.5 * (Mx - 1)) / (My - 1);
184: x[j][i].u = PetscCosReal(2.0 * PetscSqrtReal(x1 * x1 + x0 * x0));
185: }
186: }
188: DMDAVecRestoreArray(da, X, &x);
189: return 0;
190: }
192: PetscErrorCode FormDiffusionCoefficient(DM da, void *ctx, Vec X)
193: {
194: PetscInt i, j, Mx, My, xs, ys, xm, ym;
195: Coeff **x;
196: PetscReal x1, x0;
199: 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);
201: /*
202: VecSetRandom(X,NULL);
203: VecMin(X,NULL,&min);
204: */
206: DMDAVecGetArray(da, X, &x);
207: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
209: for (j = ys; j < ys + ym; j++) {
210: for (i = xs; i < xs + xm; i++) {
211: x0 = 10.0 * (i - 0.5 * (Mx - 1)) / (Mx - 1);
212: x1 = 10.0 * (j - 0.5 * (My - 1)) / (My - 1);
214: x[j][i].epsilon = 0.0;
215: x[j][i].beta = 0.05 + 0.05 * PetscSqrtReal(x0 * x0 + x1 * x1);
216: }
217: }
219: DMDAVecRestoreArray(da, X, &x);
220: return 0;
221: }
223: PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal ptime, Field **x, Field **xt, Field **f, void *ctx)
224: {
225: PetscInt i, j;
226: PetscReal hx, hy, dhx, dhy, hxdhy, hydhx, scale;
227: PetscScalar u, uxx, uyy;
228: PetscScalar ux, uy, bx, by;
229: Vec C;
230: Coeff **c;
231: DM cdm;
234: PetscObjectQuery((PetscObject)info->da, "coefficientdm", (PetscObject *)&cdm);
235: DMGetNamedLocalVector(cdm, "coefficient", &C);
236: DMDAVecGetArray(cdm, C, &c);
238: hx = 10.0 / ((PetscReal)(info->mx - 1));
239: hy = 10.0 / ((PetscReal)(info->my - 1));
241: dhx = 1. / hx;
242: dhy = 1. / hy;
244: hxdhy = hx / hy;
245: hydhx = hy / hx;
246: scale = hx * hy;
248: for (j = info->ys; j < info->ys + info->ym; j++) {
249: for (i = info->xs; i < info->xs + info->xm; i++) {
250: f[j][i].u = xt[j][i].u * scale;
252: u = x[j][i].u;
254: f[j][i].u += scale * (u * u - 1.) * u;
256: if (i == 0) f[j][i].u += (x[j][i].u - x[j][i + 1].u) * dhx;
257: else if (i == info->mx - 1) f[j][i].u += (x[j][i].u - x[j][i - 1].u) * dhx;
258: else if (j == 0) f[j][i].u += (x[j][i].u - x[j + 1][i].u) * dhy;
259: else if (j == info->my - 1) f[j][i].u += (x[j][i].u - x[j - 1][i].u) * dhy;
260: else {
261: uyy = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy;
262: uxx = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx;
264: bx = 0.5 * (c[j][i + 1].beta - c[j][i - 1].beta) * dhx;
265: by = 0.5 * (c[j + 1][i].beta - c[j - 1][i].beta) * dhy;
267: ux = 0.5 * (x[j][i + 1].u - x[j][i - 1].u) * dhx;
268: uy = 0.5 * (x[j + 1][i].u - x[j - 1][i].u) * dhy;
270: f[j][i].u += c[j][i].beta * (uxx + uyy) + scale * (bx * ux + by * uy);
271: }
272: }
273: }
274: PetscLogFlops(11. * info->ym * info->xm);
276: DMDAVecRestoreArray(cdm, C, &c);
277: DMRestoreNamedLocalVector(cdm, "coefficient", &C);
278: return 0;
279: }
281: /*TEST
283: test:
284: 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
285: nsize: 16
286: requires: !single
287: output_file: output/ex29.out
289: TEST*/