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 PetscErrorCode ierr; 118 DM da,cda; 119 120 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 121 CHKERRQ(TSCreate(PETSC_COMM_WORLD, &ts)); 122 CHKERRQ(TSSetType(ts,TSARKIMEX)); 123 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 124 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)); 125 CHKERRQ(DMSetFromOptions(da)); 126 CHKERRQ(DMSetUp(da)); 127 CHKERRQ(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 128 129 CHKERRQ(DMDASetFieldName(da,0,"u")); 130 CHKERRQ(DMCreateGlobalVector(da,&x)); 131 132 CHKERRQ(TSSetDM(ts, da)); 133 134 CHKERRQ(FormInitialGuess(da,NULL,x)); 135 CHKERRQ(DMDATSSetIFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,PetscReal,void*,void*,void*,void*))FormIFunctionLocal,NULL)); 136 137 /* set up the coefficient */ 138 139 CHKERRQ(DMDACreateCompatibleDMDA(da,2,&cda)); 140 CHKERRQ(PetscObjectCompose((PetscObject)da,"coefficientdm",(PetscObject)cda)); 141 142 CHKERRQ(DMGetNamedGlobalVector(cda,"coefficient",&c)); 143 CHKERRQ(DMGetNamedLocalVector(cda,"coefficient",&clocal)); 144 145 CHKERRQ(FormDiffusionCoefficient(cda,NULL,c)); 146 147 CHKERRQ(DMGlobalToLocalBegin(cda,c,INSERT_VALUES,clocal)); 148 CHKERRQ(DMGlobalToLocalEnd(cda,c,INSERT_VALUES,clocal)); 149 150 CHKERRQ(DMRestoreNamedLocalVector(cda,"coefficient",&clocal)); 151 CHKERRQ(DMRestoreNamedGlobalVector(cda,"coefficient",&c)); 152 153 CHKERRQ(DMCoarsenHookAdd(da,CoefficientCoarsenHook,NULL,NULL)); 154 CHKERRQ(DMSubDomainHookAdd(da,CoefficientSubDomainRestrictHook,NULL,NULL)); 155 156 CHKERRQ(TSSetMaxSteps(ts,10000)); 157 CHKERRQ(TSSetMaxTime(ts,10000.0)); 158 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 159 CHKERRQ(TSSetTimeStep(ts,0.05)); 160 CHKERRQ(TSSetSolution(ts,x)); 161 CHKERRQ(TSSetFromOptions(ts)); 162 163 CHKERRQ(TSSolve(ts,x)); 164 165 CHKERRQ(VecDestroy(&x)); 166 CHKERRQ(TSDestroy(&ts)); 167 CHKERRQ(DMDestroy(&da)); 168 CHKERRQ(DMDestroy(&cda)); 169 170 ierr = PetscFinalize(); 171 return ierr; 172 } 173 174 /* ------------------------------------------------------------------- */ 175 176 PetscErrorCode FormInitialGuess(DM da,void *ctx,Vec X) 177 { 178 PetscInt i,j,Mx,My,xs,ys,xm,ym; 179 Field **x; 180 PetscReal x0,x1; 181 182 PetscFunctionBeginUser; 183 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)); 184 185 CHKERRQ(DMDAVecGetArray(da,X,&x)); 186 CHKERRQ(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 187 188 for (j=ys; j<ys+ym; j++) { 189 for (i=xs; i<xs+xm; i++) { 190 x0 = 10.0*(i - 0.5*(Mx-1)) / (Mx-1); 191 x1 = 10.0*(j - 0.5*(Mx-1)) / (My-1); 192 x[j][i].u = PetscCosReal(2.0*PetscSqrtReal(x1*x1 + x0*x0)); 193 } 194 } 195 196 CHKERRQ(DMDAVecRestoreArray(da,X,&x)); 197 PetscFunctionReturn(0); 198 199 } 200 201 PetscErrorCode FormDiffusionCoefficient(DM da,void *ctx,Vec X) 202 { 203 PetscInt i,j,Mx,My,xs,ys,xm,ym; 204 Coeff **x; 205 PetscReal x1,x0; 206 207 PetscFunctionBeginUser; 208 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)); 209 210 /* 211 ierr = VecSetRandom(X,NULL); 212 CHKERRQ(VecMin(X,NULL,&min)); 213 */ 214 215 CHKERRQ(DMDAVecGetArray(da,X,&x)); 216 CHKERRQ(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 217 218 for (j=ys; j<ys+ym; j++) { 219 for (i=xs; i<xs+xm; i++) { 220 x0 = 10.0*(i - 0.5*(Mx-1)) / (Mx-1); 221 x1 = 10.0*(j - 0.5*(My-1)) / (My-1); 222 223 x[j][i].epsilon = 0.0; 224 x[j][i].beta = 0.05+0.05*PetscSqrtReal(x0*x0+x1*x1); 225 } 226 } 227 228 CHKERRQ(DMDAVecRestoreArray(da,X,&x)); 229 PetscFunctionReturn(0); 230 231 } 232 233 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal ptime,Field **x,Field **xt,Field **f,void *ctx) 234 { 235 PetscInt i,j; 236 PetscReal hx,hy,dhx,dhy,hxdhy,hydhx,scale; 237 PetscScalar u,uxx,uyy; 238 PetscScalar ux,uy,bx,by; 239 Vec C; 240 Coeff **c; 241 DM cdm; 242 243 PetscFunctionBeginUser; 244 CHKERRQ(PetscObjectQuery((PetscObject)info->da,"coefficientdm",(PetscObject*)&cdm)); 245 CHKERRQ(DMGetNamedLocalVector(cdm,"coefficient",&C)); 246 CHKERRQ(DMDAVecGetArray(cdm,C,&c)); 247 248 hx = 10.0/((PetscReal)(info->mx-1)); 249 hy = 10.0/((PetscReal)(info->my-1)); 250 251 dhx = 1. / hx; 252 dhy = 1. / hy; 253 254 hxdhy = hx/hy; 255 hydhx = hy/hx; 256 scale = hx*hy; 257 258 for (j=info->ys; j<info->ys+info->ym; j++) { 259 for (i=info->xs; i<info->xs+info->xm; i++) { 260 f[j][i].u = xt[j][i].u*scale; 261 262 u = x[j][i].u; 263 264 f[j][i].u += scale*(u*u - 1.)*u; 265 266 if (i == 0) f[j][i].u += (x[j][i].u - x[j][i+1].u)*dhx; 267 else if (i == info->mx-1) f[j][i].u += (x[j][i].u - x[j][i-1].u)*dhx; 268 else if (j == 0) f[j][i].u += (x[j][i].u - x[j+1][i].u)*dhy; 269 else if (j == info->my-1) f[j][i].u += (x[j][i].u - x[j-1][i].u)*dhy; 270 else { 271 uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy; 272 uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx; 273 274 bx = 0.5*(c[j][i+1].beta - c[j][i-1].beta)*dhx; 275 by = 0.5*(c[j+1][i].beta - c[j-1][i].beta)*dhy; 276 277 ux = 0.5*(x[j][i+1].u - x[j][i-1].u)*dhx; 278 uy = 0.5*(x[j+1][i].u - x[j-1][i].u)*dhy; 279 280 f[j][i].u += c[j][i].beta*(uxx + uyy) + scale*(bx*ux + by*uy); 281 } 282 } 283 } 284 CHKERRQ(PetscLogFlops(11.*info->ym*info->xm)); 285 286 CHKERRQ(DMDAVecRestoreArray(cdm,C,&c)); 287 CHKERRQ(DMRestoreNamedLocalVector(cdm,"coefficient",&C)); 288 PetscFunctionReturn(0); 289 } 290 291 /*TEST 292 293 test: 294 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 295 nsize: 16 296 requires: !single 297 output_file: output/ex29.out 298 299 TEST*/ 300