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)

106: {
107:   TS  ts;
108:   Vec x, c, clocal;
109:   DM  da, cda;

111:   PetscFunctionBeginUser;
112:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
113:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
114:   PetscCall(TSSetType(ts, TSARKIMEX));
115:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
116:   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));
117:   PetscCall(DMSetFromOptions(da));
118:   PetscCall(DMSetUp(da));
119:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));

121:   PetscCall(DMDASetFieldName(da, 0, "u"));
122:   PetscCall(DMCreateGlobalVector(da, &x));

124:   PetscCall(TSSetDM(ts, da));

126:   PetscCall(FormInitialGuess(da, NULL, x));
127:   PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, PetscReal, void *, void *, void *, void *))FormIFunctionLocal, NULL));

129:   /* set up the coefficient */

131:   PetscCall(DMDACreateCompatibleDMDA(da, 2, &cda));
132:   PetscCall(PetscObjectCompose((PetscObject)da, "coefficientdm", (PetscObject)cda));

134:   PetscCall(DMGetNamedGlobalVector(cda, "coefficient", &c));
135:   PetscCall(DMGetNamedLocalVector(cda, "coefficient", &clocal));

137:   PetscCall(FormDiffusionCoefficient(cda, NULL, c));

139:   PetscCall(DMGlobalToLocalBegin(cda, c, INSERT_VALUES, clocal));
140:   PetscCall(DMGlobalToLocalEnd(cda, c, INSERT_VALUES, clocal));

142:   PetscCall(DMRestoreNamedLocalVector(cda, "coefficient", &clocal));
143:   PetscCall(DMRestoreNamedGlobalVector(cda, "coefficient", &c));

145:   PetscCall(DMCoarsenHookAdd(da, CoefficientCoarsenHook, NULL, NULL));
146:   PetscCall(DMSubDomainHookAdd(da, CoefficientSubDomainRestrictHook, NULL, NULL));

148:   PetscCall(TSSetMaxSteps(ts, 10000));
149:   PetscCall(TSSetMaxTime(ts, 10000.0));
150:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
151:   PetscCall(TSSetTimeStep(ts, 0.05));
152:   PetscCall(TSSetSolution(ts, x));
153:   PetscCall(TSSetFromOptions(ts));

155:   PetscCall(TSSolve(ts, x));

157:   PetscCall(VecDestroy(&x));
158:   PetscCall(TSDestroy(&ts));
159:   PetscCall(DMDestroy(&da));
160:   PetscCall(DMDestroy(&cda));

162:   PetscCall(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;

174:   PetscFunctionBeginUser;
175:   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));

177:   PetscCall(DMDAVecGetArray(da, X, &x));
178:   PetscCall(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:   PetscCall(DMDAVecRestoreArray(da, X, &x));
189:   PetscFunctionReturn(PETSC_SUCCESS);
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;

198:   PetscFunctionBeginUser;
199:   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));

201:   /*
202:   ierr = VecSetRandom(X,NULL);
203:   PetscCall(VecMin(X,NULL,&min));
204:    */

206:   PetscCall(DMDAVecGetArray(da, X, &x));
207:   PetscCall(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:   PetscCall(DMDAVecRestoreArray(da, X, &x));
220:   PetscFunctionReturn(PETSC_SUCCESS);
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;

233:   PetscFunctionBeginUser;
234:   PetscCall(PetscObjectQuery((PetscObject)info->da, "coefficientdm", (PetscObject *)&cdm));
235:   PetscCall(DMGetNamedLocalVector(cdm, "coefficient", &C));
236:   PetscCall(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:   PetscCall(PetscLogFlops(11. * info->ym * info->xm));

276:   PetscCall(DMDAVecRestoreArray(cdm, C, &c));
277:   PetscCall(DMRestoreNamedLocalVector(cdm, "coefficient", &C));
278:   PetscFunctionReturn(PETSC_SUCCESS);
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*/