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