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