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