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 { 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(0); 64 } 65 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; 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(0); 102 } 103 104 int main(int argc, char **argv) 105 106 { 107 TS ts; 108 Vec x, c, clocal; 109 DM da, cda; 110 111 PetscFunctionBeginUser; 112 PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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)); 120 121 PetscCall(DMDASetFieldName(da, 0, "u")); 122 PetscCall(DMCreateGlobalVector(da, &x)); 123 124 PetscCall(TSSetDM(ts, da)); 125 126 PetscCall(FormInitialGuess(da, NULL, x)); 127 PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, PetscReal, void *, void *, void *, void *))FormIFunctionLocal, NULL)); 128 129 /* set up the coefficient */ 130 131 PetscCall(DMDACreateCompatibleDMDA(da, 2, &cda)); 132 PetscCall(PetscObjectCompose((PetscObject)da, "coefficientdm", (PetscObject)cda)); 133 134 PetscCall(DMGetNamedGlobalVector(cda, "coefficient", &c)); 135 PetscCall(DMGetNamedLocalVector(cda, "coefficient", &clocal)); 136 137 PetscCall(FormDiffusionCoefficient(cda, NULL, c)); 138 139 PetscCall(DMGlobalToLocalBegin(cda, c, INSERT_VALUES, clocal)); 140 PetscCall(DMGlobalToLocalEnd(cda, c, INSERT_VALUES, clocal)); 141 142 PetscCall(DMRestoreNamedLocalVector(cda, "coefficient", &clocal)); 143 PetscCall(DMRestoreNamedGlobalVector(cda, "coefficient", &c)); 144 145 PetscCall(DMCoarsenHookAdd(da, CoefficientCoarsenHook, NULL, NULL)); 146 PetscCall(DMSubDomainHookAdd(da, CoefficientSubDomainRestrictHook, NULL, NULL)); 147 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)); 154 155 PetscCall(TSSolve(ts, x)); 156 157 PetscCall(VecDestroy(&x)); 158 PetscCall(TSDestroy(&ts)); 159 PetscCall(DMDestroy(&da)); 160 PetscCall(DMDestroy(&cda)); 161 162 PetscCall(PetscFinalize()); 163 return 0; 164 } 165 166 /* ------------------------------------------------------------------- */ 167 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; 173 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)); 176 177 PetscCall(DMDAVecGetArray(da, X, &x)); 178 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 179 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 } 187 188 PetscCall(DMDAVecRestoreArray(da, X, &x)); 189 PetscFunctionReturn(0); 190 } 191 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; 197 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)); 200 201 /* 202 ierr = VecSetRandom(X,NULL); 203 PetscCall(VecMin(X,NULL,&min)); 204 */ 205 206 PetscCall(DMDAVecGetArray(da, X, &x)); 207 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 208 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); 213 214 x[j][i].epsilon = 0.0; 215 x[j][i].beta = 0.05 + 0.05 * PetscSqrtReal(x0 * x0 + x1 * x1); 216 } 217 } 218 219 PetscCall(DMDAVecRestoreArray(da, X, &x)); 220 PetscFunctionReturn(0); 221 } 222 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; 232 233 PetscFunctionBeginUser; 234 PetscCall(PetscObjectQuery((PetscObject)info->da, "coefficientdm", (PetscObject *)&cdm)); 235 PetscCall(DMGetNamedLocalVector(cdm, "coefficient", &C)); 236 PetscCall(DMDAVecGetArray(cdm, C, &c)); 237 238 hx = 10.0 / ((PetscReal)(info->mx - 1)); 239 hy = 10.0 / ((PetscReal)(info->my - 1)); 240 241 dhx = 1. / hx; 242 dhy = 1. / hy; 243 244 hxdhy = hx / hy; 245 hydhx = hy / hx; 246 scale = hx * hy; 247 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; 251 252 u = x[j][i].u; 253 254 f[j][i].u += scale * (u * u - 1.) * u; 255 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; 263 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; 266 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; 269 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)); 275 276 PetscCall(DMDAVecRestoreArray(cdm, C, &c)); 277 PetscCall(DMRestoreNamedLocalVector(cdm, "coefficient", &C)); 278 PetscFunctionReturn(0); 279 } 280 281 /*TEST 282 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 288 289 TEST*/ 290